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ABSTRACT 

We present 3D numerical simulations of the early evolution of long-duration gamma- 
ray bursts in the collapsar scenario. Starting from the core-collapse of a realistic pro- 
genitor model, we follow the formation and evolution of a central black hole and cen- 
trifugally balanced disk. The dense, hot accretion disk produces freely-escaping neutri- 
nos and is hydrodynamically unstable to clumping and to forming non-axisymmetric 
(m = 1, 2) modes. We show that these spiral structures, which form on dynamical 
timescales, can efficiently transfer angular momentum outward and can drive the high 
required accretion rates 0.1 — 1 M s _1 ) for producing a jet. We utilise the smoothed 
particle hydrodynamics code, Gadget-2, modified to implement relevant microphysics, 
such as cooling by neutrinos, a plausible treatment approximating the central object 
and relativistic effects. Finally, we discuss implications of this scenario as a source of 
energy to produce relativistically beamed 7-ray jets. 



Key words: gamma-rays: bursts 
neutrinos - black hole physics 



accretion disks - hydrodynamics - instabilities 



1 INTRODUCTION 



Gamma-ray bursts (GRBs) w ere first reported by 
iKlebesadel. Strong fc Olson] l| 19731 ) using data from the Vela 



satellites. The study consisted of 16 bursts lasting < 30 s, 
which were detected in an energy range of 0.2-1.5 MeV. 
Since then, the number of observed GRBs has risen into the 
thousands, using multiwavelength instruments, with burst 
durations spanning more than 5 orders of magnitude. Initial 
estimates of the total output of electromagnetic energy have 
been scaled down to ~ fo 50-52 erg with the assumption of 
relativistic beaming in jets. Observations to date of temporal 
duration and spectral hardness ratios confirm a bimodal dis- 
tribution for the population: short duration bursts (SGRBs, 
< 2 s) and long duration bursts (LGRBs, > 2 s), with dis- 
tinct progenitor scenarios assumed for each class. Here, we 
examine the latter with the aid of numerical simulations. 
In constructing realistic models, we first consider the obser- 
vational constraints on progenitors and the microphysical 
processes which control the early phases of LGRB evolution 
and provide the required conditions for jet production. 

A major breakthrough for the understanding of LGRBs 
came with their observed association with rather ener- 
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getic Type Ic (core-collapse) supernovae (SNe), some- 
times calle d hypernovae (HNe) : first, with GRB980425 and 
SN1998bw dGalama et al.lll998l), a nd then with GRB030329 
and SN2003dh jstanek et alj|2003h . Properties of these and 
subsequent LGRB-SN associations are listed in Table [T] 
Studies of host galaxies have sho wn that LGRBs ten d 
to occur in star forming regions (|Le Floc'h et all l2003l l. 
which may tend to h ave low-metallicity (< 0.3 — 0.5 Zq) 
|Fruchter et all 120061 ). The emergent picture of the LGRB 
event is that of a massive (short-lived), stellar progenitor 
collapsing in o n itself, which led to the adoption of the 'col- 
lapsar' model l|Wooslevill993l ) originally used in the context 
of core-collapse SNe. 

Briefly, the typical scenario begins with a rotating, mas- 
sive star (A//M0 > 25 — 30, possibly smaller if a remnant 
from a merger) collapsing to form a black hole (BH) from 
its core and a centrifugally-supported accretion disk, with 
continued infall from the surrounding stellar interior. The 
precise means of converting the energy of the disk mate- 
rial into high-energy photons and forming a relativistic jet 
is uncertain; however, the most studied methods, such as 
the Blandford-Znajek (B-Z) mechanism creating Poynting- 
flu x dBlandford fc Znaieklll977l : iLee. Wiiers fc Browrjlioool : 
[LI I2OOO1) and neutrino annihilation accelerating plasma 
l|Paczvnski|[l99cl ; iPopham. Wooslev fc Frve3ll999l ). require 
extreme and broadly similar conditions for the disk-BH 
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Table 1. LGRB-HNe associations, with quantities calculated from observat i ons jPian et al.l 1 20061 : 
iMazzali et"alll2006l ; iNomoto et~al"1l2006l; IValenti et aTll2008l; IChornock et al.ll2010l; IStarline et alj|2010h . 



GRB 


Distance 


Tgo 




SN 


Bol. lum. 


M( 56 Ni) 






(Zest) 


(s) 


(10 50 erg) 




(10 erg s ) 


(M©) 


(10 J km s 1 ) 


980425 


0.0085 


34.9 


0.01 


1998bw 


8.3 


0.5-1.0 


27 


030329 


0.169 


22.8 


180 


2003dh 


9.1 


0.35 


28 


031203 


0.105 


37.0 


0.26 


20031wt 


11.0 


0.55 


18 


100316D 


0.059 


~2300 


> 0.39 


2010bh 




* 


26 


060218* 


0.033 


~2100 


0.62 


2006aj 


5.8 


0.2 


26 










2003jd 


5.6 


0.3 


>15 



t first observation 2 weeks after explosion; * determined to be XRF; * values unpublished, to date. 



system in order to reproduce observed GRB energies. In 
either case, accretion rates must be quite high, Mace > 
0.1 — 1.0 Mq s _1 . Structurally, the disk must be dense 
> 10 8 " 10 g cm" 3 ) and hot (T > 10 10 K ~ 1 MeV), plac- 
ing it in a regime where various scattering processes lead 
to large amounts of neutrino production. Most regions of 
these disks, while optically thick to electromagnetic radia- 
tion, are optically thin to neutrinos and therefore are able to 
cool quite efficiently; such disks are referred to as 'neutrino- 
dominated accretion flows' (NDAFs). 

It is a great challenge for a physically realistic col- 
lapsar model to form such a disk (see discussion below 
in §2). The progenitor must possess enough specific angu- 
lar momentum so that material will balance centrifugally 
outside of the BH innermost stable circular orbit (ISCO), 
which, for a moderately spinning 2-3 Mq BH, is given by 



ied tori related to SGRB s producing high accretion and 
neutrino production rates JSetiawan. Ruffert fc Jankall2004l ; 



Jisco > 2 - 4 x 10 10 cm" 



It is then difficult to transfer 



that momentum outward efficiently to reach the required 
BH accretion rates. We note that some collapsar models 
have empl oyed a rapidly rotating magnetar as a central en- 
gine (e.g. lUsovll 19921; lThompsonf[l993 ; IWheeler et ai]|200d ; 



iBucciantini et all 20091 ). but we do not consider these alter- 
native scenarios here. 

Early studies treated NDAFs as thin disks using 
the Shakura-Sunyaev Q-viscosity ansatz to driv e the 
accretion, for either steady s t ate |Popham et al. 



iNaravan. Piran fc Kumarl l200ll; Kohri & Mincshigej 1200 



or time-dependent ( Janiuk et al.l 2004) solutions, typically 
assuming that the mag neto-rotational instability (MRI) 
l|BaIbus fc Hawlevl Il991) could provide a physical mecha- 
nis m for angular momentum transfer. A rece nt 2D study 
bv lLopez-Camara. Lee fc Ramirez-Ruij (|2009l ) with sophis- 
ticated equation of state (EoS) and cooling prescrip- 
tion found high accretion rates (> 0.1 Mq s _1 ) and 
near-GRB neutrino luminosities (< 10 52 erg s _1 ) for 
disks driven solely by the a-viscosity, though (as pre- 
sented here) the inclusion of azimuthal dynamics seems 
to be necessary to understand the full behaviour of col- 
laps ar flows. Ideal magneto- hydrodynami c (MHD) stud- 
ies l|Proga et al.1 120031 ; iMasada et al.l 120071 ) have also pro- 
duced high accretion rates, as have 2D numerical sim- 
ulatio ns studying viscosity-driven co nvection (for SGRB 
disks, [Lee. Ramirez-Ruiz fc Page|[2005l ); the creation of col- 
lapsar jets predominantly by the interplay of magnetic 
fields and neutrino annihilat ion has been invest i gated us- 
ing 2 D MHD and GRMHD |Nagataki et al.ll2007l ; iNagatakil 
120091 ). 3D simulations including general relativity have stud- 



IShibata. Sekiguchi fc Takahashill2007l ). However, the domi- 
nant physical processes during the evolution of a realistic, 
3D collapsar system remains uncertain. 

We consider the picture in which the hydrodynamic 
properties of a collapsar disk provide a mechanism for the 
high mass accretion rates giving rise to the LGRB. Since a 
BH formed from the iron core of a collapsed star has a mass 
of ~ 2 — 3 Mq, for a disk with sufficient mass to accrete 
at 0.1 Mq s -1 , the resulting disk/central object mass ratio, 
> 1/10, is already quite large compared to most astrophys- 
ical systems. This ratio, combined with the high tempera- 
tures and densities mentioned earlier, has strong implica- 
tions for the evolution of the system. Such self-gravitating 
disks are susceptible locally to clumping and are unstable 
globally to the formation of structures such as spiral arms, 
which are quite efficient at transferring angular momentum 
outwards and matter inwards. Moreover, neutrino produc- 
tion from the cooling material provides a possible mecha- 
nism for creating a jet from acceler ated plasma, in addition 
to further destabilising such disks |Rice et al.l l2003). 

It is important to note that LGRBs exhibit a wide va- 
riety of characteristics in lightcurve shapes and afterglow 
spectra , as well as in associ ations with HNe (e.g., see discus- 
sion bv lNomoto et al.ll2007l ). To date, only four LGRBs have 
been linked with HNe (Table[T]), while many partnerless ones 
have been seen at relatively close distances. Broadlined SNe 
Ic bright enough t o be classified as HN e have been observed 
without LGRBs (|Modiaz et al.1 I2008L and within); and at 
least one 'HN' has been associated with a l ess-energetic X- 
ray fl ash (XRF), XRF 060218-SN 2006aj (|Campana et al.1 
2006). Certainly, viewing effects due to jet coll imation play a 
role i n some of the differences between events |Mazzali et al.l 
2005). For example, XRFs and unaccompanied HNe may be 
explained as LGRBs viewed slightly and greatly off-angle, 
respectively (the latter interpretation is su pported by ob- 
served rates, e.g. iPodsiadlows ki et al.ll2004l ). 

It may be possible to explain much of this heterogene- 
ity by considering the variation of a few physical factors 
among (qualitatively) similar progenitors: different values 
of rotation, core size or metallicity; varied circumstellar en- 
vironments; etc. Certainly, LGRBs and related events have 
been observed across a wide range of epochs and environ- 
ments. Assuming that (most) characteristics have reason- 
able ranges for permitting bursts, one would expect to have 
not only significant variation in burst behaviour but also 
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the occurence of closely related, non-GRB events. For in- 
stance, systems with most (but not all) of the appropri- 
ate progenitor factors would likely yield 'near-' or imper- 
fect LGRBs, i.e., producing weak or unobservable jets, or 
yielding sub-critical amounts of nucleosynthetic material for 
powering HNe. Complete understanding of LGRBs would 
include both mapping out the phase space of these factors, 
as well as determining properties of related phenomena by 
altering factors of successful events. 

Here, we present a self-consistent model for powering 
LGRB engines via the collapsar mechanism using 3D nu- 
merical simulations. Starting with a physically realistic pro- 
genitor (from a stellar evolution code) soon after collapse, 
we follow the subsequent formation and evolution of an ac- 
cretion disk around a compact object (a proto-neutron star 
(PNS) which accretes matter and becomes a BH). The aim 
of this study is to include as many as possible of the impor- 
tant microphysical and dynamical features of the system. 
Therefore, it has been necessary to use several approxima- 
tions such as a simplified equation of state, appropriate for a 
radiation dominated polytrope; a time-independent estima- 
tion of electron fraction; and assumptions of nuclear statis- 
tical equilibrium (NSE), where appropriate, in determining 
constituent matter. To include effects of general relativity 
(GR) tractably, we first turn to standard pseudo-potentials 
with kinematic rescaling. We discuss the individual validity 
of these simplifications below. We note that, even while in- 
cluding several different types of physics, we generally find 
that the average estimated error for each approximation is 
< 10-20%. 

In §2, observational constraints on LGRBs are dis- 
cussed, leading to our choice of progenitor and requirements 
of disk behaviour; in §3, we briefly introduce the numer- 
ical method used here, smoothed particle hydrodynamics 
(SPH); in §4, the further physics added to the SPH code 
are outlined, including a central object, relativistic effects, 
the calculation of electron fraction and constituent nuclei in 
the material, and cooling from neutrino production; in §5, 
the initial conditions of the various models are described; in 
§6-7, results of simulations with the shellular and cylindri- 
cal angular momentum profiles, respectively, are presented, 
with those of low temperature (shellular) models in §8; in 
§9, a test scenario with a simple jet is given; in §10, com- 
parisons are made between models utilizing Newtonian and 
pseudo-Newtonian potentials; and finally, general trends and 
conclusions are discussed in §11. 



2 PROGENITOR AND ENGINE 
REQUIREMENTS 

In this section we review the broad constraints which LGRB 
observations and previously understood physics require from 
any physical collapsar progenitor. We also discuss how the 
success of any model is determined here, by judging the 
results of the formation and evolution of the disk-compact 
object system. 

2.1 Observational and physical constraints 

First, the progen itor must b e massive, M/Mq 25 — 30 
for a single-star ljFrveJll999h . in order to form a BH from 



mainly core material. Also, a large supply of matter must be 
present, assuming that rapid accretion continues for most of 
the LGRB duration and that the event may be accompanied 
by a SN; e.g., an average accretion rate of M = 0.5 Mq s _1 
for a 100 s burst would require 50 Mq, plus the mass of 
the SN ejecta as well. This progenitor range is in accor- 
dance with that of the large remnant mass (~ 3 Mq) ob- 
served for SN 1998bw and the associated SN kinetic ener- 
gies (related to th e amount of radioactiv e nickel produced) , 
2 - 5 x 10 52 erg l|Nakamura et al.1 120011 ). To date, all HNe 
associated with LGRBs have possessed broad-line spectral 
features, suggesting expan sion velocities ^ 20,000 km s _1 
(e.g., see lPian et al.l l|2006h and references within). 

Second, the LGRB-associated SNe were classified as 
Type Ic, due to the observed spectra containing neither H 
nor He lines. Therefore, it seems that, in general, the progen- 
itor should possess no envelope at the time of collapse or, at 
least, one with only a very small amount of H and He. Fur- 
thermore, this lack of an envelope assists jet propagation 
near breakout from the star, greatly reducing the amount 
of radiation energy through dissipation into kinetic energy 
(avoiding the so-called 'baryon mass-loading problem'). 

In order to form a disk, a progenitor must be rotating 
rapidly, so that the entire star does not collapse directly into 
a BH. Instead, a large fraction becomes centrifugally bal- 
anced around a PNS /BH formed by the core, which in turn 
accretes material. Disk-BH systems are common in astro- 
physics (quasars, X-ray binaries, etc.) and tend to produce 
collimated outflows aligned along or near the rotation axis. 
Such systems are also typically very efficient in converting 
the energy of infalling matter into outgoing, beamed radi- 
ation. Any centrifugally supported disk must have enough 
specific angular momentum to remain outside of the BH 
^isco, noted before as j « a few x 10 16 cm 2 s _1 . 

Finally, any collapsar disk must transport angular mo- 
mentum outwards and mass inwards with great enough ef- 
ficiency to produce the accretion rates requi red to power a 
GRB jet. In a study bv lPopham et ail l|l999T ). most models 
required accretion rates of M > 0.1 Mq s _1 to produce a jet 
by B-Z or of M > 0.1 — 1 Mq s _1 for ^-production and an- 
nihilation to produce luminosities > 10 51 erg s , although 
this depended on several factors such as BH spin, disk vis- 
cosity, etc. In the collapsar, a disk will be maintained by the 
continually infalling, high-j stellar material which reaches 
centrifugal balance. 

These few requirements are perhaps surprisingly dif- 
ficult to combine within a self-consistent, robust model. 
As mentioned previously, during the evolution of the cen- 
trifugally formed disk, it becomes a formidable challenge 
to transfer angular momentum outwards efficiently enough 
to meet the necessary accretion rates. Similarly, difficulties 
arise in satisfying the characteristics of the progenitor itself; 
a massive star can easily lose its H/He envelope, but the typ- 
ical mechanism of radiatively driven winds then ca uses the 
star t o lose angular momentum efficiently as well (|Langerl 
ll99Sl ; lLanger. Heger fc Garcfa-Seguralll99Sft . 

Here, we address the latter difficulty by utilising a pre- 
collapse progenitor evolved from the merger of two helium 
stars (more details below in fj5j) . The remnant object re- 
mains massive, while the merging process is able to eject 
outer envelope material from the system during a common- 
envelope phase l|Podsiadlowski. Joss fc Hsu 1992). The po- 
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lar direction is particularly cleared of H/He, benefitting jet 
propagation, and the remnant can easily retain sufficient an- 
gular momentum to form a centrifugally balanced disk. This 
is neither a very exotic nor overly restrictive requirement 
for the progenitor population, as at least half of all mas- 
sive s tars are estimated to be in binary and multi-star sys- 
tems l|Garmanv. Contl fe Massevl[l980l ; lKobulnickv fc Frverl 
l2007h . 



2.2 Disk evolution 

We discuss here the dynamics controlling the evolution of 
the collapsar disk in different phases. First, after the initial 
collapse, viscous processes (along with neutrino cooling) de- 
termine transport and energy dissipation. Then, given the 
large self-gravity and the rapid cooling rate of the disk, we 
expect spiral arms to form on a dynamical timescale due to 
hydrodynamic instability. While viscous dissipation is still 
present in this latter phase, the spiral structure dominates 
the global transfer of angular momentum and leads to rapid 
accretion. 



2.2.1 Disk viscosity 

The evolution of thin and slim disks has been widely stud- 
ied, typically in a low density regime but with an increasing 
interest in high density disks. Viscous processes in accretion 
disks convert kinetic energy to heat and transfer angular mo- 
mentum outwards, though generally not at rates nearly as 
high as those required to drive LGRB production. Typically, 
this viscosity is attributed to the formation of convective 
currents, due to some combination of Maxwell (magnetic) 
and Reynolds (turbulent) stresses. Of ten in studies, these 
proces ses are parametrised as the a of IShakura fc Sunvaevl 
(|l973h , similar to an efficiency, which relates the overall vis- 
cosity to the local sound speed, c s , and scale height, H: 

u = ac s H. (f) 

Though a topic of some dispute, it is reasonable to assume 
that a ^ 1, corresponding to a convection cell with size less 
than the disk height and with subsonic turnover. 

Studies of Newtonian disks have shown that the pres- 
ence of magnetic fields can lead to a magnetorotational 
instability (MRI), which ind uces efficient angular m omen- 
tum transfer and dissipation jBalbus fc Hawlevlll99ll 'l , dom- 
inating disk dynamics with a maximal growth rate on 
the order of the dynamical timescale. This analysis as- 
sumes a Boussinesq (incompressible) flow, where the rota- 
tion profile is monotonically decreasing with radius, and 
that the plasma /3 ^ 3 (ratio of hydrostatic to mag- 
netic pressure). Simulations using ideal (infinite conduc- 
tivity) MHD have shown both axisymmetric and non- 
axisymmetric perturbations, leading to large accretion rate s 
|Fromang. Balbus fc De Villiersll2004l ; iFromang et~ai1l2004 ) . 
In a physical disk, however, both the finite conductivity 
and resistivity of material can be important in damping the 
growth of the MRI, as discussed in Appendix A. We do not 
include any magnetic effects here, except by proxy through 
the influence of the a-viscosity. 

It has been shown that the numerical viscosity terms 
in SPH behave very similarly to a shearing viscosity in ro- 
tating systems. In the case of thin disks, (2D) studies have 



shown that the effective visco sity, vsph, is proportional to 
the a-viscosity (|Murravlll99rj ; Taylor & Miller, in prepara- 
tion). Therefore, the physical viscous process are essentially 
parametrised by the SPH viscosity. 



2.2.2 Hydrodynamic instability 



The Toomre parameter, Qt. (ISafronov||l960l ; lToomrdll964l ; 
iGoldreich fc Lvnden-Bel3 ll965T l has been used to charac- 
terise stability in several gaseous, stellar and hydrodynamic 
self-g ravitating disk systems (e.g.. lGammie|[200ll ; lRice et alj 
2003!) . Moreover, it has been applied to describe both local 
stability against clumping and fragmentation, and the more 
global formation of spirals and non-axisymmetric modes. 
These structures can dramatically change the evolution of a 
disk. The Toomre parameter is given by 



Qt 



ttGE ' 
1 d(a R 7? 3 ) 
R 3 dR 



(2) 



where k is the epicyclic frequency (= fl, orbital velocity, for 
Newtonian cases, and then often the stability parameter is 
approximated by Q' T = S7c s /7tGEd); art, radial acceleration; 
c s , local sound speed; G, gravitational constant; and E, sur- 
face density of the disk. The exact critical value below which 
a disk is unstable to non-axisymmetric perturbations varies 
with the specific system but is usually w 1. For example, 
a finite-thic kness, isothermal fluid disk has a critical value , 
Qt ~ 0.676 l|Goldreich fc Lvnden-Bellll965l ; lGammIe]|200ll ). 

Physically, the Toomre parameter describes a competi- 
tion between the locally stabilising influences of pressure, P 
(Cg oc T oc P/p), and epicyclic motion, and the destabilising 
influence of self gravity, which tends to enhance local dumpi- 
ness and global non-axisymmetry . From a global perspec- 
tive, lHachisu. Tohline fc Eriguchil (|l987t ) have shown that 
(for cloud fragmentation) a parameter similar to Qt can be 
derived by comparing relevant timescales in the system: the 
free-fall time, rg = (37r/32Gp) 1//2 ; the sound crossing time, 
r s d = 2-kR/cs, for a cylinder of characteristic radius, R; and 
the epicyclic period, r cp = 2n / Besides identifying two sta- 
bility criteria from the ratios Tff /r s< j and t& / r cp , the product 
of these ratios is directly proportional to the Toomre param- 
eter as given by Eq. [2] above, Qt oc Tg/T s dT cp . 

Studies of rotating polytropic disk models have sug- 
gested that m = 1 an d m = 2 instabilities can develop o n 
dynamical timescales l|lmamura. Durisen fc Pickett! [2000). 
There, angular momentum was transported outward in an 
initial pulse, followed by continued mass drain through the 
bar/spiral's interaction with outer material. 

Additionally, cooling plays a large role in the disk dy- 
namics and evolution. Besides causing vertical contraction, 
the rate of cooling, Q~, has an important effect on the 
dynamical formation of non-axisymmetric structure . Sim- 
ulations of both local and global stability l|Gammiel 1200 ll ; 
iRice et al]|2003l . respectively) showed that when the cooling 
timescale for a fluid element of specific internal energy, e, is 
of the order of the disk orbital period, i.e. 



t coo i = - — ss a few fi 1 , 

Q- 



(3) 



then massive disks can rapidly become unstable to fragmen- 
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tation and/or to non-axisymmetric structure formation on 
dynamical timescales. 

One can think of this cooling effect entering into the 
Toomre criterion via thermal pressure and the 1 st Law of 
Thermodynamics. If a fluid element in equilibrium is com- 
pressed without a cooling mechanism, then P (oc ep) in- 
creases, resisting compression and returning the volume to 
equilibrium. If the same element is compressed in the pres- 
ence of a cooling mechanism, the internal energy necessarily 
increases less (or decreases) , and P cannot restore the same 
equilibrium. In some cases, these local overdensities propa- 
gate as spirals, becoming global modes. 

While useful predictions of disk stability can be made 
from these analytic arguments, the resulting phenomena are 
highly nonlinear. Furthermore, in the collapsar case the fur- 
ther effects of relativistic dynamics (on both the epicyclic 
frequency and the disk structure) are uncertain. Therefore, 
any discussion of disk behaviour and evolution in the pres- 
ence of non-axisymmetric perturbations requires 3D numer- 
ical modelling for a range of initial conditions. 



2.3 Energetic requirements 

The full details of the process of producing a relativistically 
beamed GRB jet are uncertain, and in this study we do not 
select one of the several proposed mechanisms. Here, we dis- 
cuss the necessary conditions given by analytic arguments 
for disk-BH evolution in both the Blandford-Znajek (B-Z) 
and neutrino annihilation scenarios, which currently are the 
leading jet mechanisms and have been widely studied. The 
very complicated LGRB central engine and process of jet for- 
mation are thus characterised in terms of a few measurable, 
'physical' quantities of the evolving system, which is obvi- 
ously a great simplification and relies on many assumptions 
for each model. In the current field of LGRB studies, these 
approximations serve as the best evaluation for a collapsar 
system to produce a burst, but they can neither preclude 
new mechanisms nor definitively describe the full behaviour 
of an individual one. 

In the B-Z mechanism, rotational energy of the disk- 
BH system is converted into Poynting-flux jets along the 
rotation axis. In addition to the BH spin and mass accretion 
rates, the flux of the outflow depends on the strength of a 
central magnetic (dipole) field which reaches outward from 
the BH. For a 3 M Q BH with a 10 15 G magnetic field, the jet 
energy is given as a function of BH spin, a (= Jc/GM 2 , the 
ratio of total BH spec ific angular momentum to mass), a nd 
accretion rate, M, as l|Thorne. Price fc MacDon ald 1986): 

£ B -z~3xl0 52 a 2 (j^pr) ergs" 1 . (4) 

In the case of neutr ino annihilation, disk models studied 
bv lPopham et al.l l| 19991 ) only reached GRB- like luminosities 
of > 10 51 erg s _1 with accretion rates of w 0.1 Mq s _1 and a 
few x 0.1 Mq s _1 for a = 0.95 and 0.5, respectively. Neutrino 
production varied greatly with BH spin and disk viscosity, 
and in general, estimations to translate L v directly into jet 
energy remain quite uncertain. The efficiency of neutrino an- 
nihilation in plasma is unknown, particularly in the presence 
of other physical factors such as magnetic fields, baryon pol- 
lution, etc. Analysing idealised disks around Kerr BHs using 



ray-tracing techniques, iBirkl et al.l (|2007f ) estimated the ef- 
fects of both GR and disk geometry on neutrino annihilation 
rates. In general the effects were small (most less than a fac- 
tor of 2), with more efficient annihilation associated with 
thinner disks and higher BH spins. For this study, we note 
in particular the important influence of the latter, as the 
smaller i?isco with increased a leads to greater L„ in the 
inner regions as well. 

The neutrino case can b e generalised more broadly with 
the simple assumption fe.g. Jjanka et al" I I1999T ) that the use- 
able neutrino energy rate for producing a LGRB jet, E v , 
scales directly with the total disk L v : 

E v = r\ v pL v , (5) 

where r) v „ now represents the total, unknown efficiency for 
the physical process of neutrino annihilation. This is a very 
uncertain quantity, and studies have shown that this effi- 
ciency may not be constant. Here, we will consi der it to have 
a reas onably hig h value, r\ V y = 0.01, as found bv lJanka et al.l 
(1999) (see also. iRuffert fc Jan ka 1998), accepting that this 
is a very rough approximation. 

In this study .Egrb = 10 51 erg s _1 = 1 foe s _1 is de- 
fined as our fiducial LGRB-jet energy flux, the threshold 
value which a successful central engine must produce. This 
value is generally consistent with observed LGRB energet- 
ics and would then require, for a BH with a = 0.5, for ex- 
ample, an accretion rate of M > 0.13 Mq s _1 in the case 
of the B-Z mechanism. However, we do not intend this to 
preclude the possibility of lower-luminosity bursts which, 
as of this time, have not been de tected (see discussion in 
IChapman. Priddev fc Tanvir| [2009). In fact, such objects or 
'failed' LGRBs may produce other observed phenomena, 
such as polarised SNe, XRFs, etc., particularly in the pres- 
ence of a large central B-field. 



3 SPH 

For thi s study we use smoothed particle hyd rodynamics 
(SPH) (|Gingold fc Monaghanl [l977l ; lLucvlll977D . which has 
found wide usage in modeling various astrophysical phenom- 
ena. SPH has advantages over other numerical techniques in 
that it naturally adapts to following dynamic flows and ar- 
bitrary geometries without the need for mesh refinement or 
other readjustment techniques required in Eulerian-based 
codes. The Lagrangian method allows natural definitions 
for the energetics of fluid elements and also efficiently maps 
structures whose density spans several orders of magnitude. 
We utilise the public release version of the Gadget-2 code 
(Springcl 2005), which, furthermore, strictly conserves mo- 
mentum (linear and angular), energy and entropy. We add 
additional relevant physics to the code, such as cooling, rel- 
ativistic effects and the presence of a central object. 

Briefly, in SPH a continuous fluid is sampled at a fi- 
nite number of points, and discretised versions of the three 
hydrodynamic equations (mass, momentum and energy) are 
calculated. Gadget-2 utilises a polytropic EoS and, rather 
than energy, evolves an entropic function, K(s) = P/p 1 , 
where 7 is the usual adiabatic index; P and p, the hydrody- 
namic pressure and density, respectively, related to specific 
internal energy, e = (7 — l)P/p; and s is specific entropy. 
The smoothed kernels determine a set of nearest neigh- 
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bours which interact with symmetric hydrodynamic forces. 
Entropy is created from discontinuities and shocks, via a 
viscosity formulation which possesses two terms. The first 
provides the bulk and shear viscosity (linear in converg- 
ing velocity), and the second, an artificial von Neumann- 
Richtmyer (quadratic) viscosity, mimics naturally dissipa- 
tive pro cesses in additio n to preventing interparticle pene- 
tration l|Monaghanlll997l ). The strength of the Gadget-2 vis- 
cosity is parametrised by a single coefficient of order unity 
(here, a v = 0.8), and the Balsara prescription to reduce nu- 
merical viscosi t y introduced in rotating flows was utilized 
ijBalsaral 1 19951 : fsteinmetj 1 19961 ). Density is not explicitly 
evolved but instead is recalculated at each timestep from a 
weighted average of the total mass within a smoothed kernel. 
Gravitational forces are calculated efficiently with a hierar- 
chic al multipole expansio n, using standard tree algorithms 
(e.g. IBarnes fe Hutlll986l ). 



4 ADDITIONAL PHYSICS 

To these general hydrodynamic and gravitational compo- 
nents, we include physical features specifically relevant to 
the coUapsar evolution. Unless otherwise stated, all temper- 
atures, T, are given in units of K, and both mass and proton 
density, p and Y c p, respectively, in units of g cm -3 . 



material, and we do not expect a strong model-dependence 
from using these relations. 

A fluid particle accretes onto the central object in two 
ways. Most basically, when it reaches Rpns with specific an- 
gular momentum, j, then the particle's mass, m, and total 
angular momentum about the axis of rotation, J = mj, are 
added to that of the PNS (calculating j is further discussed 
in £|4. 2|) . Also, since SPH particles represent a density dis- 
tributed through a kernel, we 'accrete' a fraction of the mass, 
mf, and total angular momentum, J = mtj, as the parti- 
cle approaches -Rpns to within a fraction of the smoothing 
length. This action also requires that any accreting particle's 
smoothing length must be much less than -Rpns, effectively 
establishing a resolution requirement for the simulation. 

Above its maximum mass, the PNS becomes a BH, 
whose properties are determined, in this case, solely by its 
mass and spin. During this stage, the flow boundary pos- 
sesses no surface pressure and is approximated as the sur- 
face of a sphere with radius n n = Risco of a Kerr BH in 
the equatorial plane. All matter is co-rotating with the BH , 
and r- m can be calculated exactly (jNovikov fc Frolovlll989l . 
and see Eq. ((7}, below). The numerical procedure of the ac- 
cretion is the same as in the PNS case, and we note that the 
final Rpns is fortuitously similar to n n for a newly-formed 
BH, which prevents the introduction of numerical errors at 
the conversion. 



4.1 Central object 

After collapse, the iron core of the initial progenitor first 
forms a (rotating) PNS, which subsequently accretes mate- 
rial to become a BH. Since our main concern is studying the 
behaviour of the accretion flow, we focus on simulating only 
the main effects of each distinct phase on the hydro-flow. 
The two most important influences of the central object are 
to provide an inner accretion boundary for the inflowing 
material and to exert gravitational accelerations (see H4.2I 
below, regarding relativity). Here, we discuss how the PNS 
and BH differ importantly in the presence and absence, re- 
spectively, of a surface pressure. 

In the first few seconds after formation, the central ob- 
ject is quite hot and is best described as a PNS, with slightly 
larger radius than the typical 10-12 km of a NS, as well 
as a 'softer' surface. As the PNS accretes matter, its ra- 
dius decreases, and we approximate its ma ss-radius (M-R ) 
progression using the recent calculations of iNico tra (2006). 
Examples (using baryonic mass) of (Mpns/Mq, .Rpns /km) 
from his model T30Y04 are: (1.6, 35), (1.7, 24) and (1.9, 
16), beyond which a BH is formed. 

To approximate the surface of a (spherical) PNS, we 
include a radially-outward acceleration at the given -Rpns, 
which interacts with SPH particles that have approached to 
within a fraction of their smoothing-length. For inflowing 
material, we add into the particle's Euler equation an accel- 
eration term which mimics the steep pressure gradient of the 
outer PNS, aR = — VP/p, approximated from the Tolman- 
Oppenheimer-Volkoff (TOV) equation for our model. How- 
ever, fluid elements with density greater than that at neu- 
tron drip (p > pdrip = 4x 10 11 g cm -3 ) accrete directly. The 
important aspect in these specifics is to provide a reason- 
able representation of the PNS effects on the evolving disk 



4.2 General and special relativity 

In order to include effects of general relativity in an ap- 
proximate manner, we add the acceleration of a pseudo- 
Newtonian potential for the central object to the gravity 
calculated between the SPH particles, in both the BH and 
PNS phases; though, the effects are larger for the BH case 
due to its larger mass and greater spin, and also Rpns > r- m 
for all masses, so that the innermost flow is at a larger ra- 
dius, further decreasing PNS-GR effects. In this section we 
utilise relativistic units with c = G = 1, and, while the same 
treatment of acceleration is used for both types of compact 
object, below we refer to them collectively as 'BH'. 

Due to the large angular momentum of the accreted 
material, one must utilise approximations appropriate for 
a rotating central object. From the many pseudo-potenials 
which have been developed in different studies, we have cho- 
sen to use the form of inward radial free-fall acceleration, 
aR,(R), for a Kerr BH with mass, M, and specif i c an gular 
momentum, a, given by iMukhopadhvav fc Misral |2003l ): 



fa\ M 
asEp(R) = -jg 



(6) 



derived from what they call the 'Second-order Expansion 
Potential' (SEP), where all radii have been scaled to grav- 
itational units, r g = GM/c , and n n is the location of the 
ISCO, given by: 



3 + Z2-[(3-Zi)(3 + Zi+2Z 2 



,1/2 



(7) 



Zi = l + (l-a 2 ) 1/3 [(l + a) 1/3 + (l-a) 1/3 ], 
Z 2 = (3a 2 + Z 2 ) 1/2 . 

We discuss this choice of pseudo-potential and compare its 
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characteristics to others, as well as to full GR solutions, in 
Appendix B. 

In order to keep track of the BH spin parameter, a, 
we add to it the total angular momentum of all accreted 
material, J = mj. For material accreting from the equatorial 
plane, we approximate the specific angular momentum, j, as 
that of a corotating particle in circular orbit at rj n , given as 
jShapiro fc Teukolsky|[l983l ): 



Jk = 



2aVMr in + a z 



3Mr in + 2a/Mf~) 1 /2 ' 



(8) 



which is a several times larger than the Newtonian equiva- 
lent in the vicinity of the accretion boundary. In the case of 
polar accretion (defined as when the distance to the equa- 
torial plane is greater than that to the rotation axis), we 
simply use the Newtonian value, j = xv y — yv x . Necessarily, 
this material will tend to have much less angular momentum 
anyway, so that this is a reasonable approximation. 

For all velocities discussed here, we utilise a simple 
special-relativistic rescaling of the values occuring within 
the pseudo-potential, v ps : 



VSR 



(9) 



This has been shown by lAbramowicz et all (|l99rj ) to pro- 
duce velocities within 5% of relativistic calculations around 
non-rotating BHs when used in conjunction wi th the 
Paczynski-Wiita potential (|Paczvriski fc Wiita| [l980). 

We recognize that using these approximations to the 
Kerr metric in order to represent the gravitational contri- 
bution of the central object is a rough way of proceeding, 
but we believe that it constitutes a legitimate approach at 
the present stage of the work (see further discussion in Ap- 
pendix B). In the case where the central object is a BH, 
it is clearly very approximate to treat the gravitational ac- 
celeration on fluid elements as the sum of contributions of a 
vacuum BH pseudo-potential and the Newtonian self-gravity 
of surrounding matter. In the case where the central object 
is a PNS, the surrounding spacetime differs from that of 
the Kerr metric anyway, with the former having quadrupole 
moments considerably larger than those for equivalent BHs. 
(However, it must be noted that the effect of the quadrupole 
falls off as R 3 , so the influence of this difference for the disk 
flow may not be very great.) Improving on this treatment 
remains a topic for future work, but we believe that the 
approach being used here is a reasonable one at present. 

4.3 Neutrino cooling 

In the inner regions of the collapsing star, the flow is op- 
tically thick to electromagnetic radiation, so that radiative 
transfer is very inefficient and the energy carried away by 
photons is negligible. Instead, much of the material, partic- 
ularly in the disk itself, occupies a region of (Y c p, T) phase 
space where cooling via neutrino emission is important (cal- 
culating V c is discussed in the next section). In the collapsar 
system local neutrino production greatly influences disk sta- 
bility, and, summed globally, it may provide a mechanism 
for relativistic jet production. 

For the regions which are optically thin to neutrinos, 
energy from processes such as pair annihilation, plasmon 
decay and (non-degenerate) bremsstrahlung is carried away 
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Neutrino Cooling Rates log(erg cm" 3 s' 1 ) 




log(g cm"') 



Figure 1. Solid lines show contours of the total combined cooling 
rate, Q v . The dark shaded regions show where a single produc- 
tion process dominates > 90% of the total cooling rate, and the 
lighter, > 75%. Included for reference is the electron degeneracy 
temperature, T Cit j cg . 



freely. Very dense regions with p > 5 x 10 12 are considered 
to be optically thick to neutrinos, and using the common 
'lightbulb approximation', cooling shuts off. 

For the ranges 7 log(T) sC 12 and 6 ^ log(Y c) o) «S 
12, we utilise tables of data from detailed calculations 
by lHaft Raffelt fc Weiss! (I 19941 ) for the plasma processes 
and by lltoh et al l j 19961 ) for pair annihilation; for non- 
degenerate bremsstrahl ung scattering, we use the analytic 
approximation given bv lHannestad fc Raffeltl (|l998l ). Q Br = 
1.5 x W 33 T?i 5 pl 3 erg cm" 3 s" 1 , with the notation T u = 
T/10 11 , etc. Fig. Q] shows contours of the combined rates of 
neutrino cooling, Q v , with shaded regions where each pro- 
cess dominates total cooling. During the evolution of the 
collapsar disk, we note that nearly all material for which 
cooling is relevant occupies (Y e p, T) states dominated by 
plasma processes and pair production. 

The importance of the relations between various 
timescales for the stability of the collapsar disk has already 
been introduced in £12.2,21 (and tnse is discussed in i )4.5l) . 
Fig. [2] shows the relevant cooling timescales (t C ooi = t/Qv) 
for neutrino processes in this phase space. Typical values of 
the dynamical timescale (td yn = fl , orbital period) around 
a 3 M@, a — 0.5 BH at different radii have also been calcu- 
lated. For direct comparison, the dynamical timescales are 
overlaid on the phase space, with reference to the cooling- 
stability criterion of Eq. <(3j . Fluid at a given radius is stable 
against clumping for the various (Y e p, T) states in the re- 
gion below the td yn (r) contour, and unstable for those a few 
times above. 
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Figure 2. A comparison of timescales: the black (white) dotted 
lines show contours of t coo i (tnse)j from Eq. (3). In most regions 



TNSE <C t 

c 



Values of dynamical timescale, 4<j y 



given and overplotted on the t coo i contours for direct comparison 
regarding local stability. 



Figure 3. Contours show the average nucleon number, A, in NSE 
for a given temperature and density, as calculated by F. Timmes 
(private communication) . The dashed curve shows the commonly 
used onset of free nucleon matter. Vertical dot-dashed lines show 
density-dependent electron fraction, Y e . 



4.4 Electron fraction 

As the density of matter increases, so does the rate of cap- 
ture of electrons by nuclei. This neutronisation decreases the 
electron fraction, Y e , and necessarily the proton/neutron ra- 
tio, below 0.5. This change has important consequences for 
nuclear processes and cooling via neutrino production, which 
are dependent on the proton density, Y c p. For densities 
p > 1 7 , we follow the simple prescription of iLiebendorfed 
|2005l ). who noted in core-collapse simulations that Y c re- 
mains nearly constant in time and temperature for a given 
density. Using their 'G15' set of parameters, the approxi- 
mate electron fraction, Y e (p) is given by: 

Y e (x) = 0.354 - 0.1 Ux 

+ 0.035|x| [1 + 4 (|^ -0.5)(|x| - 1)] , 

x(p) = max [-l,min(l, 0.3434 log(p) - 3.5677)] . (10) 

The electron fraction for a range of densities is shown by the 
vertical lines in Fig. [3] Thus, roughly-speaking, the electron 
fraction remains both temperature- and time-independent, 
with the restriction that, for a given ith fluid element, Y Bt i 
monotonically decreases over time, i.e., upwards fluctuations 
are not permitted. Matter with density p < 10 7 retains its 
initial value of Y e ,i. 

4.5 NSE 

In regions of high temperature and density, nuclear reac- 
tions involving strong interactions equilibrate with their in- 
verse reactions, and nuclear statisti cal equilibrium (NSE) 
is reached (e.g.. lilix fe Meyer] I20061 . and within). Assum- 
ing that matter remains transparent to neutrinos, ele- 
ment/isotope abundances may then be calculated using 



statistical relations which involve a large network of rel- 
evant isotopes, for material with T ^ Tnse (= 4 x 
10 9 K » 0.34 MeV). Importantly, in these regions Fe- 
group elements are formed, predominantly from Si-group 
elements created durin g a lower (p,T), quasi-static equilib- 
rium phase (QSE. e.g. |Bodansky. Clayton fc Fowlen 1 19681 : 
iMever. Krishnan fc Clavtonlll998l ) 

For calculating the average nucleon number, Anse, of 
a fluid element during the simulation, we consider matter 
with both T > T NSE and p > 10 6 to be in NSE (for a 
discussion of tnse, see below). Computationally, a fluid el- 
ement below Tnse retains its pre-collapse A value. Know- 
ing (p, Y c , T), tInse is determined from a table of values 
calculated by F. Timmes (private communication, and see 
lTimmes| [l999) and shown as a function of temperature and 
density (as above, Y c — Y c (p)) in Fig. [3] Also shown as a 
guide is the often-used approximation for where the mass 
fraction o f free neutrons and pro tons in material is unity, 
X nuc = 1 (|Wooslev fc Baronlll992l ). 

We calculate the abundances of specific a- and heavy- 
elements produced by NSE conditions in the inner regions 
of the collapsar using results of the same reaction network 
calculations by F. Timmes (see Appendix C). Due to the 
observational association of LGRBs with HNe/SNe, we in- 
vestigate in particular the production of 56 Ni, the radioac- 
tive isotope which is the dominant energy source of SN 
lightcurves via decay to 56 Co. A typical core collapse SN 
contains 56 Ni(~ 0.1) at time of explosion, and estimates for 
the brighter HN may require up to 0.5—0. 7 Mq of 56 Ni; some 
SNe m ay possibly require up to ~ 5 M© l|Umeda fc Nomotq 
(2008) and references within), but this latter value remains 
uncertain. The masses of competing end products, non- 
radioactive 54 Fe and elemental Co, are also approximated; 
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free nuclei and light a-elements are expected to dominate 
much of the hot, dense collapsar disk. 

Elemental abundances provide further means for un- 
derstanding the observed association of asymm etric and po- 
larised SNe with LGRBs. liviazzali et all (|2005h showed that, 
among other features, a HN viewed along the axis of the 
GRB jet would possess a strong, single-peaked [Fe II] emis- 
sion line, while one seen in the equatorial plane of the system 
would possess a double-peaked [O I] emission line. It was 
reasoned that Fe-group elements were produced during the 
collapse near the jet, while the lighter elements, remaining 
from the progenitor's evolution, settled into the equatorial 
plane and were ejected from the evolving disk. Their analysis 
matched the Fe-line observed in SN 1998bw (associated with 
GRB 980425), and the O- feature in SN 2003jd (no GRB). 

Therefore, we examine the amounts and locations of Fe 
and O produced in the collapsar from NSE, noting the dif- 
ference in origin for the latter in this case; at larger radii, 
a significant mass of O created during the stellar evolu- 
tion of the progenitor is expected to remain in roughly a 
spherical distribution. The timescale for reactions to reach 
NSE varies strongly with temperature and very weakly 
with density. Here, we utilise the approximation tnse = 
exp(149.7/T 9 -46.054) s (|Gamezo. Khokhlov fe Oranll2005h . 
This value for fluid elements is compared with other rele- 
vant timescales in the disk (shown in Fig. particularly 
the dynamical and viscous timescales in the inner and outer 
parts, respectively. In regions where the nuclear timescale is 
significantly less than either of these, one can calculate the 
nucl ear abundances of m aterial with relative ease (again, see 
e.g., lHix fe MeverH2006T ). 



4.6 The role of non-NSE nucleosynthesis 

The 56 Ni required to power HN lightcurves may also be pro- 
duced by non-NSE nucleosynthetic tracks, for example dur- 
ing the expansion and cooling of hot material. Such time- 
dependent nuclear processes (e.g., the r-process and explo- 
sive nucleosynthesis) have not been included in this study. 
However, the feasibility of producing a significant amount 
of 56 Ni by non-NSE processes can be gauged through an 
examination of the properties of the material comprising 
candidate regions for rapid expansion, namely in outflows 

(i) from the outer edge of disk: equatorial outflows due to 
angular momentum transfer; 

(ii) from equatorial plane: vertical neutrino or thermally 
driven winds from hot disk material; 

(iii) from the inner edge of the accretion disk: hot 'bub- 
bles' coming from non-accreted flow or created from jet- 
interaction, ejected in the polar direction. 

In typical cases for SNe and outflows from compact ob- 
jects, Fe-group and heavier elements are manufactured by 
high entropy material (of order 10 he, (A and higher, as a 
minimum, where fee is Boltzmann's constant) which is ex- 
panding rapidly (v ~ 0.1c). Furthermore, the electron frac- 
tion of material (Y e < 0.5 initially for most regions of in- 
terest here) is important in determining nucleosynthetic re- 
sults; outflow with an extremely low electron fraction per- 
mits r-processing (yielding neutron-rich elements), while 
that which is more in the neighbourhood of 0.44 < Y c < 0.5 



is a candidate to evolve to Y c w 0.5, at which point sig- 
nificant 56 Ni may be produced (though the precise evo- 
lution depends on multi ple timescales, on the presence of 
neutrinos, etc.; see, e .g. Pruet. Wooslev fe Hoffmanl 2003 : 
Jordan fe Meved |2004 iPruet. Thompson fe Hoffmanl I2OO4I ; 
Surman. McLaughlin fe H ix 2006; Scit enzahl et al.ll2008l ). 



5 PROGENITOR MODEL 

The 'physical' progenitor used in these simulations comes 
from the resultant of the merger of 2 He st ars (8+8 M©), 
evolve d to the point of core collapse by iFrver fe Hegerl 
| 2005j ). We utilise the results of their model 'M88c', which 
was predicted to form a BH instead of a NS and which also 
possessed a large amount of angular momentum (the max- 
imum of their several scenarios of the merger process). We 
now desecribe the mapping of the ID, pre-collapse (poly- 
tropic, 7 = 5/3) model into our 3D axially symmetric, post- 
collapse (7 = 4/3) starting point, and also the range of 
related models which we analyse here. 

The collapsar simulation begins ~ 2 s after core col- 
lapse, after the stellar core has formed a 1.75 Mq PNS 
(> Mchandra = 1.4 M© , but at this stage, the PNS is quite 
'oversized' before neutrinos have diffused out). A rarefac- 
tion wave moving radially outward at the local sound-speed 
leads to co ntinued infall. Followin g Colgate. Herant fe Bend 
l|l993h and lFrver. Benz fe HerantJ j 19961 ). we assume that an 
isentropic, convective envelope in pressure equilibrium has 
formed around the PNS. This material possesses a high spe- 
cific entropy per baryon (s cnv = 15 Its /A), and at its outer 
radius (w 10 7 cm) the hydrostatic pressure balances the ram 
pressure (=pv 2 ) of the infalling stellar matter. 

From the lD lFrver fe Hegerl (120051 ) profile, we first gen- 
erate a 3D pre-collapse model which is spherically symmetric 
in density and other EoS variables. Similarly, we assume that 
Ye and A (not necessarily in NSE) are spherically symmetric. 
We implement two general classes of rotation profile: first, 
shellular, where the angular velocity, (j>, is constant for a 
given radius from the centre; and second, cylindrical, which 
is expected for a state of permanent rotation with aligned 
surfaces of constant density and pr essure, as coro llaries of 
the Poincare-Wavre Theorem (e.g.. iTassoull Il978l ). Shellu- 
lar rotation is expected in systems with strong, anisotropic 
turbulence oriented perpendicular to equipotential surfaces 
|Zahnlll992l ). 

From the above, there are three distinct regions to be 
translated from the pre-collapse progenitor: PNS, envelope 
and infall. The inner 1.75 Mq simply becomes the PNS, 
whose outer radius and numerical properties have been de- 
scribed in ij4.ll The envelope is then built up in spherically 
symmetric layers of hydrostatic equilibrium, with radial den- 
sity profile: 



P (R) 



GMpt 



RpiSIS 



1/4 



(11) 



where Ppns is the pressure at the inner edge of the envelope. 
The entropic function in the envelope, K cnv = P/p 1 , defines 
the thermal profile and is, for a radiation-dominated gas, 
simply related to specific entropy per baryon, s, as K = 
(3c/4<r) 1//3 (fcB s/4m p ) 4//3 , where a is the Stefan-Boltzmann 
constant, and m p , the proton mass. 
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Table 2. Summary of collapsar models, mapped from IPrver fc Hegerl {2005) pre-collapse binary merger (details in text). Rotation 
and E t h profiles have been scaled in cases to test a wide array of scenarios. Values of quantities are given for the initial Mo = 0.5 Mq 
of free material in the models, which utilise the SEP pseudo-potential unless otherwise stated. At t = s, the total potential energy 
for the SEP and Newtonian potentials is E^j ggp = —11.9 X 10 50 erg and -Ev.Ncw = —12.0 X 10 50 erg, respectively; all models have 
the same total radial kinetic energy, -Ekin R = 4.3 X 10 50 erg, and the integrated (inward) radial mass-flux, Mr = 3.44 Mq s — 1 . 



Model 


Velocity 


<p 




Jo 


_ 


-^kin [-^kin,</> ] 




Further 




Profile 


Scaling 


Scaling 


(10 50 erg s) 


(jo / jiQpn! 


(10 50 erg) 


(10 50 erg) 


Notes 


PSm2 


shell 


x2 


full 


15.64 


65.5 


7.78 [3.48] 


3.83 


- 


PSml 


shell 


full 


full 


7.82 


32.8 


6.04 [1.74] 


3.83 


- 


PSdsq2 


shell 


/V2 


full 


5.53 


23.1 


5.17 [0.87] 


3.83 




PSdsq2J* 
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^The azimuthal contribution to -Ekm- ^Exactly the same evolution as PSdsq2 until t = 1.10 s. 



To determine the structure of the infalling material, we 
use an approximate rate of fallback, M < 1 Mq s _1 , across 
a spherical shell a t 2 x 10 7 cm , appropriate for collapsing 
massive stars (e.g- lFrveifeoOrj ). The radial velocity is taken 
to be that of freefall, vg, from the original radius of the shell, 
and therefore, the infalling density is given by pi n t(R) = 
M/(4ivR 2 Vf{). For the azimuthal velocity, v$, we preserve the 
specific angular momentum of the fluid elements in collapse. 
Finally, the internal energy of the fluid elements is calculated 
by assuming nearly steady flow with negligible cooling, so 
that the Bernoulli equation applies. 

In both the envelope and the infalling region, for mate- 
rial with p > 10 6 and T > 10 7 (T > 4 x 10 9 ), Y c (A NS e) is 
calculated explicitly, as described in H4.4I Otherwise, those 
quantities are unchanged in a fluid element from their pre- 
collapse values. The PNS (and later, BH) surface defines 
the inner boundary condition. At the outer boundary mass 
is regularly added to the system as continued infall, ensuring 
that the outermost radius at any time is much larger than 
the accretion disk/inner region of interest. At t = s, the 
simulation contains Mo = 0.5 M@ of free material modelled 
using 2.25 x 10 SPH particles of equal mass (which remains 
the same for all added infall particles). 

Table [2] lists the properties of models tested in this 
study, calculated only from the initially present Mo. They 
do not include properties of the subsequent infalling layers, 
and as such, the bulk values are relevant measures of relative 
energies and angular momentum, but they are not absolute 
quantities for the entire collapsar. In addition, we define the 
dimensionless ratio, yj, of the inflow average specific angu- 
lar momentum, jo = Jo /Mo, to the characteristic quantity, 
jisco of the central object (even for the Newtonian poten- 
tials, as an approximation), and we discuss further applica- 
tion of this in £1111 In order to further explore a wider range 
of possible LGRB progenitors, various physical properties of 
the pre-collapse progenitor are rescaled. In particular, we fo- 
cus on quantities which appear in Qt, affecting the stability 
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R/UO 8 cml 

Figure 4. Model PSml: evolution of the Toomrc parameter, Qt- 
The disk remains globally stable, as Qt > 1- One local ring 
(where Qt ~ 1) develops non-axisymmetry which eventually is 
damped globally (see Fig. l5t. 



of the forming accretion disk. Scalings of both the shellu- 
lar and cylindrical rotation profiles are tested, which affect 
E by determining the location of centrifugal balance for a 
fluid element. Similarly, the specific internal energy (oc c a ) 
is decreased, to examine the effects on M and the neutrino 
cooling rates by changing the initial location of a fluid ele- 
ment in the (p, T) phase space. For comparison, a Newtonian 
potential is used in two models to gauge the effects of rela- 
tivity on the collapsar system. Finally, in one model a 'jet' 
is inserted at late times after the formation of a spiral, in 
order to examine the production of outflows. 
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PSml: log(p), a=0.438, M BH =1.8071, t=0.30s PSml: log(p), a=0.457, M BH =1.8355, t=1.00s 
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Figure 5. Model PSml (equatorial density) closely follows the evolution of Qt-, as the disk is globally stable to hydrodynamic perturba- 
tions present in the early phases. A locally unstable region leads to the growth of temporary non-axisymmetric structure, which damps 
out due to global stability 



6 SHELLULAR ROTATION MODELS 

We now describe the evolutions of various models with shel- 
lular profiles. In order to appreciate the evolution of a num- 
ber of the more dynamic models, we also make reference to 
animation j]] showing the evolution of density (p), specific 
entropy per baryon (s/A) and entropic function (K). 

6.1 PSml 

Fig. [4] shows the evolution of Qt, for the evolving PSml 
accretion disk. Also shown are estimates of the mass of the 
forming disk compared to central object mass (Md/Mq). As 
the disk grows, Qt approaches the region of instability in 

1 Located at: | htt p : / / www- astro . physics . ox . ac . uk/ ~ pt aylor | All 
animations show the evolution of a thin slice of the equatorial 
plane, excepting PSdsq2_invden_mov.gif (vertical plane slice). 



one localised region. At late times, the minimum remains 
just greater than unity throughout the disk, which is there- 
fore predicted to remain globally stable against the devel- 
opment of non-axisymmetric structures (though locally, one 
region remains near instability). In testing the often-used 
approximation for classifying stability, it was found that the 
shape of the Q' T curve was very similar to that in Fig. [4] 
but with systematically smaller values, in roughly constant 
ratio, Q't/Qt « 70 — 80%. This behaviour remains in the 
other pseudo-potential systems as well, and only Qt has 
been used in analysis. 

The disk behaviour during evolution, as shown by the 
equatorial density contours in Fig. [5J closely follows the 
Toomre stability analysis. Large non-axisymmetric pertur- 
bations which appear at early times (t ~ 0.30 s) have be- 
come globally axisymmetric by t = 1.00 s. The small, m = 2 
mode in density which is beginning to form at R < 10 7 cm 
(at t = 1.00 s), where Qt ~ 1 predicted a locally unstable 
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PSml: log(T), t=1.00s 



PSml: Y„ t=1.00s 




Figure 6. Model PSml properties at t = 1.00 s: T and Y e in the equatorial plane (top); in the vertical plane, p and s (middle), and Y c 
and A (bottom), for which in the latter are included lines of A = 1 (dotted), 2 (dashed), 3 (dot-dashed), 4 (dot-dot-dashed). 



© 2008 RAS, MNRAS 000.ITH291 



Disk instabilities in LGRBs 13 



FSml: v R (cm s" 1 ), t=1.00s 



PSml: Accretion Rate (M a s"') 



10'. 



10 ; 



10 s 



1.0 1.5 2.0 

R/(10 8 cm) 



Figure 7. Radial velocity profile of PSml at t = 1.00 s, for which 
'diamonds' and '+' show outflows and inflows, respectively. 



region, develops into a higher order mode by t = 1.50 s. 
While outer regions become temporarily disturbed, the disk 
remains globally stable and returns to a quasi-steady stat^l 
by t = 2.00 s. Thus, the instability remains local and does 
not lead to the formation of large spirals and to enduring, 
non-axisymmetric modes. The disk is dense (p > 10 8-12 ) 
and hot (T > a few MeV), with inner regions reaching elec- 
tron fractions Y c < 0.3, as shown in the top panel of Fig. [6] 

The vertical disk structure is shown in the lower panels 
of Fig. [(J] The density and specific entropy show the inner 
'thin disk' region extending out to Rt < 2 x 10 7 cm, with an 
extended quasi-thin disk out to « 2 Rt- A large 'corona' of 
hot (T ~ 1 MeV) and fairly dense material extending from 
the equatorial plane is in vertical near-equilibrium (i.e., not 
in an outward- driven wind, but instead, with relatively small 
velocities) and balances the continued infall, with uniform 
Y e and dominated by material with A ^ 4. At irregular in- 
tervals, radial outflows are produced in the equatorial plane 
(Fig. [3 bottom), with velocities up to ~ 1000 km s - , which 
carry the low-A,-Y e material outwards. 

The rate of accretion, M, onto the central object is 
shown in Fig. [5] along with the disk neutrino luminosity, 
L„. During early times, M is large (while low angular mo- 
mentum material accretes quasi-spherically) , as is L v since 
infalling material shocks around the PNS surface and be- 
gins to cool rapidly, with a maximum just as M begins to 
decline. By late times, a stable disk has formed with steady 
(predominantly equatorial) accretion, and neutrino produc- 
tion decreases. By Eq. Q, these values of M are at least a 
factor of 3-4 too low to produce Eb-z ~ £grb- 

For neutrinos, even utilising a reasonably optimistic an- 



2 The disk is not an isolated system; the collapsing star provides 
continued infall (at a slowly decreasing rate), and the properties 
of successive shells vary according to the pre-collapse progeni- 
tor structure. Therefore, 'quasi-steady' refers to an axisymmctric 
flow in which profiles of density, velocity and other hydrodynamic 
variables change slowly with time. 




Disk Neutrino Luminosity (foe s" 1 ) 




Figure 8. Model PSml: evolution of mass accretion rates and 
neutrino luminosity. After a brief period of rapid accretion of 
low-j material, M decreases rapidly as the disk reaches a steady, 
axisymmctric state. The neutrino luminosity behaves similarly; 
after an initial period of rapid cooling of dense, shock-heated ma- 
terial, L„ decreases below 1 foe s _1 . 



PSml: Relative abundances, O/tO+Fe), t=1.00s 




-0.5 0.0 0.5 

x/(10 8 cm) 



Figure 9. For model PSml, contours compare the relative loca- 
tions of O and Fe produced by the disk and inner region in NSE 
(vertical slice above the equatorial plane) . In the central (empty) 
region, amounts of both elements are negligible. 



nihilation efficiency, r) vii L v <C -Egrb- Therefore, this model 
will not produce a successful LGRB by either mechanism, as 
accretion rates from the usual viscous processes by which the 
thin a-disk evolves are too low. The collapsing material pos- 
sesses too much angular momentum and thermal support to 
be globally unstable, even with large amounts of neutrinos 
being produced from the inner regions. 

To analyse the collapsar from the point of view of pro- 
ducing a HN, we have calculated the quantities of various 
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elements produced in the disk and inner regions in NSE, 
which are given for PSml (and for other models) in Table [3] 
As noted above, light elements dominate the disk material, 
and most of the heavy metals are produced in the corona and 
shocked infall regions. For PSml, the NSE mass of 56 Ni re- 
mains fairly constant at 0.02 M©, which is nearly double the 
mass of 54 Fe (0.012 Mq) and roughly equal to total elemen- 
tal Fe (0.021 Mq, though the iron content slowly increases 
with time). Therefore, NSE-nucleosynthesis in and around 
the (steady) collapsar disk falls short of the amount of 56 Ni 
required to power the bright HNe by roughly an order of 
magnitud^E 

A comparison of the relative production of elemental Fe 
and O in the fluid elements is shown in Fig. [9] As the disk 
settles into a steady flow, there seems to be only a slight pref- 
erence for Fe being produced in the polar direction and O 
equatorially in the coronae around the dense disk. It should 
be noted that what is plotted is only the relative, local abun- 
dances, and that the total mass of O produced in this case is 
several orders of magnitude less than the Fe. However, this 
mass difference may be unsurprising given that Fe, unlike O, 
is an NSE end product which dominates the mass-fraction 
for a large area of relevant (p, T) phase space (cf. Appendix 
C). 

We now investigate the three regions of potential non- 
NSE nucleosynthesis as mentioned in c]4.6l For region (i), 
some equatorial outflow material reached velocities of order 
10 3 km s _1 at R ~ 10 s cm, where Y c « 0.43 - 0.45 (Fig. ij). 
The entropy through much of the disk is quite low, but in 
the outermost regions and corona reaches ~ 10 ks/A. How- 
ever, efficient 56 Ni-production would probably require signif- 
icantly greater ejection velocities as well as higher entropy. 

There is not a great deal of outflow vertically from the 
disk for region (ii). The coronal material is optically thin to 
neutrinos and possesses a quasi-steady vertical structure; the 
velocities (not shown) of fluid elements are relatively small 
(<C 0.1 c) and non-uniform, except in the azimuthal direc- 
tion (for near-Keplerian rotation). Therefore, these coronal 
regions do not appear to be a possible site for r-process nu- 
cleosynthesis, and these characteristics are similar through- 
out the disks analysed here. 

Finally, for region (iii) near the inner edge of the disk, 
some material has Y c w 0.3, which is beneficial for re- 
processing, though the specific entropy in these regions is 
quite low (s/A ~ 1) due to efficient neutrino cooling. How- 
ever, as seen in Fig. [6] high-s regions surround the inner disk 
edge, where material is still less neutron rich Y c w 0.43—0.45. 
Furthermore, unlike in region (ii) where vertical motion of 
material is inhibited by continued infall, after the initial 
core-collapse the cone around the rotation axis has a rel- 
atively low density. One might expect vertical trajectories 
in this region either from material shock-heated to form a 
high entropy 'bubble' which rises, or from material swept 
up in the collimated outflow which creates the actual LGRB 
(discussed briefly in 
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Figure 10. Model PSdsq2: evolution of the Toomre parameter, 
Qt- By t = 0.70 s, a significant fraction of the disk is unstable 
(Qt < 1): an d spiral structure develops dynamically. 

6.2 PSdsq2 

We now investigate a model {v$ globally decreased by a fac- 
tor of \/2) with one half of the -Ekin,^ of the marginally unsta- 
ble PSml. The evolution of Qt is shown in Fig. llQI as well as 
the disk/central object mass ratios, which are significantly 
larger than in the previous model. A large part of the inner 
disk has Qt < 1 and is therefore predicted to be unstable to 
density perturbations. Indeed, non-axisymmetric structures 
form by t = 0.50 s and Fig. [TT] (top) shows spirals in the 
density evolution as well as localised density clumps. The 
same figure shows the effect of shocks in transferring kinetic 
energy into thermal, easily seen in the T-profiles. The (on- 
line) animations of density and specific entropy illustrate the 
full dynamics of the global behaviour of the disk, in particu- 
lar the creation of oscillating, m = flows in the equatorial 
plane. 

The mass accretion and neutrino production rates are 
shown in Fig. 1121 with the same initial phase of low-j ac- 
cretion. Even after the initial onset of spirals at t > 0.50 s, 
M remains low until t > 1.10 s, at which time spikes of 
rapid accretion develop. In this case accretion rates with 
M > 0.6 M Q s _1 give rise to Eb-z « 5 foe s -1 > Egbb- L„ 
increases by an order of magnitude during the rapid accre- 
tion to « 30 foe s _1 . Therefore, while the system provides 
conditions for a central engine of sufficient power to create a 
GRB jet via the B-Z mechanism, the process of neutrino an- 
nihilation is probably too inefficient to produce a 'successful' 
jet, as E u ~ 0.3 foe s _1 < -Egrb- 

The size of the T > 1 MeV inner region of this collap- 
sar is larger than that of PSml, and a much smaller mass of 



3 However, this simulation encompasses only a single accretion 
event. It is likely that continued infall would lead to disk refor- 
mation, and subsequent repetitions of accretion events may then 
itcratively increase the total 56 Ni yield. 



4 Due to the combined effects of the high accretion rates, den- 
sities and rapidly changing dynamics, the timesteps of the SPH 
particles decreased rapidly in these late times, causing the simula- 
tion to halt. These numerical considerations are discussed further 
in [ITT] 
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PSdsq2: log(p), a=0.469, M BH =1.8564, t=0.70s PSdsq2: log(p), a=0.508, M BH =1.8915, t=1.25s 
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PSdsq2: log(T), t=0.70s PSdsq2: log(T), t=1.25s 
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Figure 11. Model PSdsq2 in the equatorial plane, non-axisymmetric evolution of p and T at t = 0.7 s (top) and t = 1.25 s (middle); 
spiral structure develops and leads to shocks and angular momentum transfer, (bottom) Electron fraction and specific entropy in the 
vertical plane at t = 1.25 s. 
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Table 3. Representative NSE-nucleosynthetic yields (in Mq) at late times for each model. The largest fluctuations in 
elemental abundance (> order of magnitude) occured in PSdsq2, for which the values are given at two times. 
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Figure 12. Model PSdsq2: mass accretion rates and neutrino 
luminosity. 



PSdsq2: k B /A, Velocity field, t=1.20s 




Figure 13. Specific entropy contours with normalised velocity 
vector field showing direction of flow in vertical slice half-plane 
(above equatorial plane) during accretion phase. The flow directly 
above the poles shows small turbulent eddies. 



Fe-group elements is present, particularly by t = 1.25 (Ta- 
ble[3|. As a consequence, only a very small amount of 56 Ni is 
created, and there is no preferential production of Fe in the 
polar direction at late times. The region of A < 4 material 
is less flattened along the equatorial plane, and likewise the 
distribution of O vs Fe is less polarised (particularly due to 
the extremely low iron content at late times). 

The coronal structure (Fig. 1111 bottom panel) is smaller 
than that in PSml. It remains fairly steady, even dur- 
ing periods of outflow, and the non-aligned nature of the 
(non-azimuthal) components of velocity vectors is shown 
in Fig. 1131 ruling out case (ii) outflows for nucleosynthesis. 
Peak equatorial outflows reach a factor of a few higher than 
for PSml, but are still probably too low for case (i). The 
low-entropy region has a smaller vertical height from the 
equatorial plane, and material with higher s/A and mod- 
erate Y e lies quite close to the central object in this model 
(Fig. lll|l . Minimum values of the electron fraction in the disk 



have decreased due to high T and p of the shocked regions 
to roughly Y c — 0.25, but the specific entropy remains low. 
Again, case (iii) outflows in the presence of a jet flow remain 
the most likely location for significant 56 Ni-production. 

6.3 PSd2 

The model PSd2 (v$ decreased by a factor of 2) possesses 
25% of the -Ekin,^> of PSml and continues the trend of 
showing greater instability at earlier times after collapse 
(t w 0.20-0.25 s). Fig.nHshows the evolution of Q T and the 
increasing disk/central object mass ratio. Fig. 1 151 shows the 
perturbations evolving as a complex configuration: tightly 
wound spirals develop into filaments within a global m = 2 
mode (full evolution in the online animations). The inner re- 
gion is comprised of alternating dense 'fingers' and shocked 
flow with higher entropy. Though intially a smaller disk 
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Figure 14. By / = 0.25 s. part of the disk has reached an unsta- 
ble state (Qt ~ 1), and a larger extent of the disk has become 
unstable by t = 0.30 s. 
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Figure 16. Model PSd2: mass accretion rates and neutrino lu- 
minosity. 



forms than in previous models, the non-axisymmetric struc- 
ture extends to larger radii. 

For much of the disk evolution, the mass accretion rate 
is consistently high fFig. I16[). After the influx of \ow-j ma- 
terial, non-axisymmetric structures form (t > 0.20 s) and 
disk accretion increases to M — 0.4 M© s - , peaking at 
the high rate of M « 1.2 M s . This is reflected also in 
the neutrino luminosity, which reaches L v > 100 foe s . 
Therefore, in this model a successful central engine for a 
GRB jet is produced by both mechanisms being examined 
in this study. High accretion rates (with a ~ 0.5 and assum- 
ing a large, central B-field) lead to Eb-z ~ Eghb, while 
the large neutrino production from the shock-heated disk 



leads to E v > -Egrb (assuming a reasonable ^-annihilation 
efficiency) . 

As in the previous models, mainly light elements con- 
stitute the disk. Although the inner regions here produce a 
much greater amount of heavy elements, 0.048 Mq of 54 Fe, 
0.069 M of Fe, 0.029 M of Co and 0.037 M of 56 Ni, 
the latter remains below that required to power typical HNe 
(though, as noted previously, a repetition of accretion events 
may additively increase abundances). The outflow velocities 
in the equatorial region reach 5 x 10 s < dr < 10 9 cm s _1 . 
Much of the material at the outer edge of the disk has high 
specific entropy with Y c ~ 0.45, making case (i) outflows 
a possible site for non-NSE nucleosynthesis. While there is 
some vertical flow from the equatorial plane, most material 
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remains confined to the coronal regions, limiting the possi- 
bility of case (ii) processes. 

There is very little asymmetry observed in Fe and O, 
with production of the latter slightly increased along the 
equatorial plane near the coronal structures. 

6.4 PSd5, PSm2 

Neither model PSd5 nor PSm2 (v^, decreased by a factor of 
5 and increased by a factor of 2, respectively) results in a 
successful LGRB progenitor. In model PSd5, the accretion 
onto the central object remains quasi-spherical, with 45% 
of M coming from the polar direction even at t = 0.20 s. 
The infall rate is quite steady, with no outflows. This model 
has too little rotation to form a significant disk structure or 
a potential LGRB, and it shows very little polarisation in 
element abundances. 

In PSm2, the Toomre parameter remains 2 at all 
times. Since PSml, with 25% of the kinetic energy of this 
model, was only locally unstable, it is unsurprising that 
this disk appears to be globally stable. The mass of dense 
(p > 10 8 g cm -3 ) disk material is significantly less than that 
of model PSml, and by t = 1 s the central object has ac- 
creted approximately half as much mass. In this case, the 
distribution of elemental O is significantly polarised along 
the equatorial plane. 

Results for both models are given in Tables O and 2] 



7 CYLINDRICAL ROTATION MODELS 

We briefly discuss the cylindrical rotation profile models, 
which have lower i?kin,0 but higher Jo than the shellular 
models. In general, the results for respective scalings of ve- 
locity for each series are qualitatively similar. 

7.1 PCml 

The evolution of Qt and of the system's disk/central object 
mass ratios are shown in Fig. [IT] for PCml (full v$). In 
this model the collapsing material forms a disk without any 
locally unstable region, even though the disk average Qt > 
10 is less than in PSml. 

The densities in the inner disk region remain approx- 
imately an order of magnitude less than in PSml, with 
a slowly decreasing density profile in the outer regions 
(Fig. I18|) . The temperature and radial velocity profiles are 
similar, although Y c remains above 0.35. In general, the ver- 
tical structure and the relative locations of Fe and O are 
similar to those of PSml. Roughly a factor of five less 56 Ni 
is produced in PCml, which produces ~ 60% more Fe and 
54 Fe instead (Table [3]). 

After a brief accretion period of low-j material, M set- 
tles to « 0.03 Mq s -1 . Similarly, even with 100% neutrino 
annihilation efficiency, E v remains below Egkb • Therefore, 
this model, like PSml, does not appear to be a candidate 
for either a LGRB or a HN. 

7.2 PCdsq2 

The evolution of Qt for PCdsq2 (v$ decreased by a factor 
of v2) is shown in Fig. 1191 The disk becomes unstable by 
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Figure 17. For PCml, Qt 2> 1 for the entire evolution, and the 
disk remains stable. 



PCml: log(p), a=0.438, M BH =1.7971, t=1.00s 
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Figure 18. Equatorial density profile for PCml, which shows 
peak density less than PSml. 

t — 0.35 s, earlier than for PSdsq2 and with a larger radial 
extent. The resulting formation of spirals and an m = 2 
structure is shown in Fig. 1201 High velocity radial outflows 
form, vb. ~ 10 9 cm s _1 , carrying material with Y c w 0.43 
outward. The values of A increase much more rapidly with 
radius than in PSdsq2, with much of the material in the 
outflow regions having A > 4. 

In the vertical plane the coronal regions extend to 
greater height, with a much larger amount of material with 
Y c « 0.43. Again, the average nucleon number increases and 
the temperature drops more quickly outside of the main 
disk/coronal area. The outer corona again contains O pref- 
erentially, with the Fe abundance dominating at larger radii. 
The amount of 56 Ni (0.096 M ) remains below HN values, 
however, by a factor of 3-5. This is similar to PSdsq2, where 
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Figure 19. Qt evolution for PCdsq2, which is unstable by t 
0.35 s. 



O was more prevalent in all regions, and significantly smaller 
amounts of Fe-group elements were produced. 

After the initial infall of low-j material, accretion rates 
(shown in Fig. 1211) remain low, M < 0.05 M© s _1 , even 
after the formation of spiral structure, although there is an 
upturn at late times. Neutrino luminosity is fairly high at 
L v « 10 foe s _1 , but E v remains below -Bgrb for reasonable 
annihilation efficiencies. 



7.3 PCd2 

As in the shellular case, PCd2 decreased by a factor of 
2) becomes unstable to the formation of non-axisymmetric 
structure quickly (t < 0.20 s, Fig. [221 ■ Fig. 1251 shows the 
spiral modes which create shocks in the fluid flow. The evo- 
lution of structure, with the development of finer spirals fol- 
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Figure 21. Model PCdsq2: evolution of mass accretion rates and 
neutrino luminosity. 



lowed by a global m w 2 mode, is very similar to that of 
PSd2. Eventually, moderate- Y a material moves outward in 
the equatorial plane due to the transfer of angular momen- 
tum in the disk. 

In the vertical plane again p, s/A and Y c decrease more 
slowly with distance from the equatorial plane than in the 
equivalent shellular model, PSd2. 56 Ni and Fe are produced 
in very similar quantities to those of PSd2. Typically, the 
outer corona contains O preferentially, with Fe abundances 
being greater at larger radii. 

During the spiral structure phase, accretion rates begin 
near M < 0.1 M s~\ increasing to M < 0.4 M s" 1 at late 
times. The neutrino luminosity is fairly high during the spi- 
ral structure phase and increases to L u — 100— 150 foe s _1 . 
This model therefore provides sufficient fuel for a LGRB cen- 



© 2008 RAS, MNRAS OOO.fTlMl 



20 P. A. Taylor, J. C. Miller, Ph. Podsiadlowski 



+ tf 



PCd2: Q T 



AOA 



0.1 1 

0.0 



+ t = 0.05, M D /M C = 0.024/1.774 = 0.014 

x t = 0.10, M D /M C = 0.062/1.783 = 0.035 

O t = 0.20, M D /M C = 0.177/1.791 = 0.099 

A t = 0.30, M n /M r = 0.325/1.799 = 0.181 



0.4 0.6 
R/(10 ! cm) 



0.1 

0.0 



PSmlKcB: Q T 





1 




, | , , , 


, , , , J- , _ 

o 


" □ 








is : 

® □ 












X 

- £ 

■ o <s 
□ a 








_% 


+ 


I = 


0.10, M D /M C = 


0.019/1.765 = 0.010 






t = 


0.30, M D /M C = 


0.349/1.821 = 0.192 







t = 


0.70, M D /M C = 


0.254/1.836 = 0.139 "_ 




A 


t = 


1.00, M D /M C = 


0.227/1.843 = 0.123 




□ 

l 


t = 


1.40, M D /M C = 

1 T T 1 


0.267/1.853 = 0.144 



0.4 0.6 
R/(10 ! cm) 



Figure 22. Qt evolution for PCd2, which is unstable by t > 
0.20 s. 



Figure 24. Qt evolution for PSmlKd5, which remains stable 
globally. 
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PSmlKd5: log(p), a=0.459, M BH =1.8427, t=1.00s 
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Figure 23. Density structure in the equatorial plane for model 
PCd2 at t = 0.44 s. 



tral engine for both the B-Z mechanism and ^-annihilation, 
as did PSd2. 



8 LOW TEMPERATURE MODELS 

The low temperature (shellular) models produced results 
very similar to those with similar rotation profiles discussed 
above in Sj6] Comparisons between the relevant models are 
discussed here. 



8.1 PSmlKd5 

The mass of the disk which forms in PSmlKd5 (full v$) 
is greater that of PSml. The shapes of the Qt profiles are 
similar, with PSmlKd5 tending to lower values, though its 




Figure 25. Density in the equatorial plane for PSmlKd5, which 
has a very similar profile to that of PSml, excepting the local 
instability (cf. Fig. [5}. 



minimum value increases after t w 0.70 s, and the disk re- 
mains globally stable throughout its evolution (Figs. |2425[) . 
Much smaller regions of outflow are produced and with lower 
Mr, due to the greater stability of the disk. 

Through the first 1.00 s, PSml and PSmlKd5 have sim- 
ilar temperature profiles in the equatorial plane. In the lat- 
ter, density decreases and specific entropy increases more 
slowly with radius, and, while the electron fraction is greater 
in the inner region (Y e w 0.30), Y e increases slowly so that 
a larger region contains material with Y c « 0.43 — 0.45; A 
behaves similarly. 

Vertically, in PSmlKd5 density decreases rapidly with 
distance from the equatorial plane, and a smaller corona 
surrounds the disk than in PSml. In this region, there is a 
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Figure 26. For PSd2Kd5, Q T becomes unstable after / = 0.20 s, 
slightly earlier than PSd2; the shapes of the curves are very sim- 
ilar, with PSd2Kd5 being « 10% less (cf. Fig. [Tit. 
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Figure 27. Equatorial density profile for PSd2Kd5, which shows 
very similar late time structure to PSd2 (cf. Fig. 1151 . 



larger extent of low s/A ^ 1 material, before the entropy 
increases rapidly. The electron fraction behaves similarly in 
the corona, with Y c mainly w 0.43. The amount of 56 Ni pro- 
duced is smaller by an order of magnitude, as are the total 
masses of the Fe-group elements calculated here in general 
(Table fJJ). The relative regional distributions of O and Fe 
are similar to those of PSml, with a slightly smaller coronal 
extent of O. 

Both M and L v closely resemble those for PSml. After 
the initial accretion of shock-heated, low-j material, accre- 
tion rates settle to M ~ 0.03 M s" 1 , and neutrino luminos- 
ity goes to L v w 0.1 foe s _1 . Therefore, model PSmlKd5 is 
also unable to produce a successful LGRB. 



8.2 PSd2Kd5 

The shapes of the Qt curves for PSd2Kd5 (v^ decreased 
by a factor of 2) shown in Fig. [26] are very similar to those 
of PSd2, with the former having » 10% lower values. For 
t < 0.20 s, the disk masses are nearly identical as well. Again, 
tightly wound spirals develop into 'finger-like' substructure 
in a global m — 2 mode (Fig. I27[) . It is perhaps surpris- 
ing that the thermal and specific entropy structures are also 
nearly identical to those of PSd2, despite an initially large 
difference in internal energy. Also, the A and Y e values are 
very similar for both models. However, the outflow veloci- 
ties in the equatorial plane are much greater in PSd2Kd5, 
reaching ur > 0.1c. 

The vertical profiles of PSd2Kd5 closely resemble those 
of PSd2, but the coronal regions are smaller in the former, 
and p decreases more rapidly with distance from the equa- 
torial plane. The relative locations of Fe and O production 
follow those in PSd2. However, the amounts of Fe-group ele- 
ments produced are quite different (Table Here, 0.10 Mq 
of 56 Ni is made, approximately three times the amount in 
PSd2 and only a factor of a few below that which is es- 
timated for some HNe. The masses of the competing end- 
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Figure 28. Model PSd2Kd5: evolution of mass accretion rates 
and neutrino luminosity. 



products, 0.017 M of 54 Fe and 0.029 M of Fe, are a factor 
of three lower than in PSd2. 

In this model accretion reaches very high rates, M « 
1 Mq s _1 at peak, with an average of M ~ 0.25 — 
0.30 Mq s -1 otherwise (Fig. 1280 , The neutrino flux reaches 
the high values of L„ > 200 foe s - , although it remains 
below 50 foe s _1 for much of the evolution. Therefore, while 
at peak values both Eb-z and E v are > Sgrb, the B-Z 
mechanism appears to be significantly better for providing 
a longer duration central engine. 
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PSdsq2J: k B /A, Velocity field, t=1.20s 



NSml: Q T 




Figure 29. Specific entropy contours with normalised velocity 
vector field showing direction of flow in vertical slice half-plane 
(above equatorial plane) for comparison with Fig. 1131 Material 
in coronal regions retains small, turbulent velocities. A strong 
vertical outflow is created by the jet. 



9 SIMPLE JET TEST, PSDSQ2J 

In order to test the behaviour of the collapsar flow in the 
presence of a collimated outflow, we implemented a very sim- 
ple 'jet' structure originating from the central object (start- 
ing with model PSdsq2). This was not intended to be a phys- 
ically realistic jet, but merely a toy model to gain insight into 
affected regions of the system, i.e. in disk structure, in the 
nature of any coronal material swept outward, etc. 

The prescription inserted into PSdsq2 at t = 1.10 s was 
that a narrow jet region with polar angular radius of = 
10° around the rotation axis developed and grew quickly. 
An acceleration directed along the envelope of the jet, away 
from the equatorial plane, was applied to every fluid element 
which approached this region to within a small fraction of 
its smoothing length. 

In Fig. [291 model PSdsq2J is shown 0.10 s after the 
insertion of the jet. A comparison with the same model 
without a jet in Fig. [13] shows that the vertical structures 
are affected only slightly, but most importantly, high s/A 
material (with moderate Y c ) is carried upwards from the 
corona and inner edges of the disk at high velocities, 0.05- 
0.1c. These conditions make the inner disk region a strong 
candidate for producing large amounts of 56 Ni in outflows. 
We note, too, that the vertical outflow material becomes 
more likely to form Fe-group elements, assuming NSE is 
maintained; therefore, the polar direction becomes more Fe- 
rich, while the equatorial plane remains O-rich, as in po- 
larised HNe. Explosive nucleosynthesis in the jet itself may 
also be a source of heavy elements such as 56 Ni, though the 
abundances which can be prod uced remain uncertain (e.g., 
iNagataki. Mizuta fc Satoll2006t ). 
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Figure 30. Evolution of Qt for the Newtonian, shellular model, 
NSml. 



10 COMPARING GR AND NEWTONIAN 
POTENTIALS 

A simple test to gauge the effects of including GR was made 
by utilising a purely Newtonian central potential in (shellu- 
lar) simulations with full rotation (model NSml) and with 
decreased by a factor of 2 (model NSd2). In general, 
the results are quite similar to those using the SEP pseudo- 
potential (PSml and PSd2, respectively). 

Fig. 1301 shows Qt for NSml. The greatest difference oc- 
curs in the minimum-region which nears instability, which is 
~ 50% wider in the pseudo-potential system. The result is 
that the NSml disk remains both globally and locally sta- 
ble, without the formation of non-axisymmetric structures 
which developed in the PSml model. Due to this (tempo- 
rary) lack of perturbation, the radial outflow is much lower 
in the Newtonian case. 

The profiles of hydrodynamic and microphysical quanti- 
ties, in both equatorial and polar directions, are very similar. 
However, NSml produces orders of magnitude less 56 Ni and 
a factor of five less Fe than PSml. The M and L v values are 
typically « 15 - 20% greater in NSml. 

For the reduced velocity case, the Qt curves for NSd2 
(Fig. 13 1 p closely resemble those of its pseudo-potential ana- 
logue, PSd2. Moreover, the structural, hydrodynamic and 
nuclear constituent properties of the evolving disk appear 
similar, as well. In this case, the only difference in Fe-group 
production appears in 56 Ni, which is 25% less than the 
amount produced by PSd2. Again, the M and L„ curves 
for NSml are qualitatively similar to those of PSml but 
with the values of each of the former being ~ 20% greater. 



11 DISCUSSION AND CONCLUSIONS 
11.1 Numerical timesteps 



It was noted in H6.2l that the timesteps in the more dynamic 
simulations became too short to continue. In Gadget-2 stan- 
dard hydrodynamical criteria assign the length of time be- 
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Table 4. Results for the various collapsar models, in terms of maximum accretion rates and energy production. The values of Eb— z.max 
have been calculated using Eq. Jit . In general, scenarios of M max > 0.1 Mq s _1 and L„ jnl!a > 100 foe s -1 produce successful LGRBs 
(although the profiles for each exhibit rapid fluctuations). 



Model Vel. E th y, M max #B-Z,max ii/.max M( 56 Ni)t Instability 

Profile scale scale (jo/jisco) (Mq s — 1 ) (foe s — 1 ) (foe s — 1 ) (Mq) (description) 



PSm2 


shell 


x2 


full 


65.5 


0.023 


0.13 


0.073 


10- 4 


none 


PCml 


cyl. 


full 


full 


36.3 


0.035 


0.20 


0.23 


0.005 


none 


PSml 


shell 


full 


full 


32.8 


0.032 


0.21 


0.41 


0.020 


local perturbation (temp.) 


PSmlKd5 


shell 


full 


/5 


32.8 


0.027 


0.17 


0.12 


0.001 


none (slight non-axis.) 


NSml 


shell 


full 


full 


32.8 


0.035 


0.241 


0.20 


io- 42 


none 


PCdsq2 


cyl. 


/V2 


full 


25.7 


0.067 


0.38 


17.5 


0.010 


global spirals, late m = 2 


PSdsq2 


shell 


IV2 


full 


23.1 


0.67 


8.33 


29.7 


10 -8.3 


global spirals, late m = 


PCd2 


cyl. 


/2 


full 


18.2 


0.37 


2.51 


167 


0.032 


global spirals and m = 1, 2 


PSd2 


shell 


/2 


full 


16.4 


1.23 


8.82 


121 


0.038 


global spirals and m = 2 


PSd2Kd5 


shell 


/2 


/5 


16.4 


1.00 


6.49 


214 


0.101 


global spirals and m = 2 


NSd2 


shell 


/2 


full 


16.4 


0.63 


4.41 


121 


0.029 


global spirals and m = 2 


PSd5 


shell 


/5 


full 


6.5 










none (no disk) 



t Estimated at late times in the evolution. 
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Figure 31. Evolution of Qx for the Newtonian, shellular model 
with reduced rotation, NSd2. 



tween 'kicks' for each SPH particle: the Courant-Friedrich- 
Levy (CFL) condition, which requires small timesteps for 
dense material (At B tc P oc p _1 ); and a dynamic condition, 
which requires small timesteps for particles with rapidly 
changing velocity (A£ s t cp oc |a| -1 / 2 ). 

In the case of collapsars, densities often reach high val- 
ues, much greater than neutron drip. Similarly, dense mat- 
ter in the inner disk is shocked by spiral arms, and infalling 
matter is shocked to high T at the edge of the corona. Fi- 
nally, particles are removed from the simulation when they 
are accreted, and during periods of rapid accretion, the com- 
puter must keep searching for and recalculating neighbours 
for nearby un-accreted particles; the material close to the 
accretion boundary is also typically the densest, meaning 
that when small amounts accrete, neighbour lists must be 
recalculated for many particles. For some simulations in this 
study, the conjunction of these conditions caused timesteps 



to become quite small, and the computational expense per 
timestep, quite large. 

In order to prevent some of these difficulties, par- 
ticularly due to neighbour searches when accreting, one 
may utilise a method of replacing many SPH particles 
in a high-density region w ith a single particle, as in 
iBate. Bonnell fc Price! (119951 ). A simplified version of this 
procedure was tested in some models with only limited suc- 
cess. To avoid these difficulties in future work, more refined 
numerical treatments of this process and of the accretion 
will be required. 



11.2 Hydrodynamic Instability 

In all of the models, the global stability of the collapsar disks 
consistently depended on the scaling of azimuthal velocity in 
the progenitor. The dimensionless parameter, yj — jo/jisco, 
was introduced to characterise the progenitor in terms of rel- 
evant quantities; here, systems with values of, very roughly, 
10 — 15 < j/j < 30 — 32 (with the upperbound seen more 
directly) were hydrodynamically unstable; for large yj, an- 
gular momentum stabilises the disk, and for small values, no 
significant disk forms. The initial velocity profile shape and 
even the scaling of Eth by half an order of magnitude had 
comparatively very little influence. Results from the simula- 
tions are summarised in Tablef3](for reference, the scalings of 
the various model progenitors are included as well). We note 
that in future studies, the quantity yj may be particularly 
useful in predicting the behaviour of a collapsar-candidate 
progenitor. 

The full rotation (and Vj, x 2) models proved to be glob- 
ally stable against forming non-axisymmetric structures, 
with PSml temporarily developing a locally unstable ring. 
By scaling v$ in the progenitor by 1/ y/2, the resulting disk 
became unstable to thin spirals which became tightly wound 
and shocked surrounding fluid. These spirals eventually led 
to the development of larger, radial oscillations (m = 
modes) for PSdsq2 and to global non-axisymmetric struc- 
ture for PCdsq2. In all cases, scaling v$ by 1/2 led to the 
rapid formation of tightly wound spirals and also to the ap- 
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Figure 32. A relative comparison of specific angular momentum 
profiles for the shellular models PSmi, PSdsq2 and PSd2 at t = 
i.00, 1.25, 0.42 s, respectively. Analogous behaviour occured in the 
cylindrical and thermally-scaled models. All values are scaled in 
arbitrary units to the Newtonian, Keplerian profile (line). 



pearance of global m = 1 or 2 modes which dominated the 
subsequent disk evolution. 

As a general trend, instability set in a number of dynam- 
ical timesteps earlier for lower -Bkin,^, with cylindrical mod- 
els (prefixed with 'PC') tending to become unstable slightly 
more quickly than the shellular ones. The spirals typically 
formed in the disk within the radial interval of 4 — 8 x 10 6 cm 
(of order 20 x GM/c 2 ), with lower rotation models becom- 
ing unstable at the smaller end of the range. We note that 
for several models, disks formed localised density clumps 
as well, with signficant repercussions for the accretion rate 
profile and energetics. 

We note that it is not surprising that the profile 
of each of the models which formed global structures and 
non-axisymmetric modes had been scaled downward from 
that of the original progenitor. The total angular momen- 
tum of a binary remnant depends heavily on the redistri- 
bution mechanism of the merger process, which remains an 
area of active research. For the basic progenitor used in this 
study, the merger prescription had maximised the possible 
J retained by a remn ant, one of three pictures presented by 
iFrver fc Hegerl (|2005l ). Therefore, it should not be worrying 
for the general consistency of binary mergers as collapsar 
progenitors that only the models with scaled-down -Bkin,^ 
produced globally unstable disks. 



11.3 Accretion rates, v luminosity 

Accretion rates were nearly constant (after accreting low-ji 
material) for the full-rotation models which did not form 
spirals, settling to be « 0.03 Mq s _1 , even in the model 
with local instability. As shown in Table [4] these rates are 
too low by a factor of 5-8 to produce a burst via the B-Z 



mechanism, and the accompanying E v is also several orders 
of magnitude below the LGRB limit. Moreover, these mod- 
els would provide no central engine mechanism for the large 
variations in observed lightcurves which occur on the or- 
der of ms timescales (though other phases of evolution may 
contribute). 

In the v$/y2 models, accretion rates only increased 
slightly after the initial onset of the spiral phase and be- 
fore the formation of global modes; in the PSdsq2 model, 
the accretion rate then briefly peaked at M > 0.65 Mq s~ . 
The neutrino luminosity was ~ 10 foe s _1 during the early 
spiral phase, increasing by a factor of two in PCdsq2 and 
by a factor of three in PSdsq2 as global modes appeared. 
While some of the accretion rates reach LGRB energetic 
requirements, the E v typically remained too low with the 
annihilation efficiency taken into account. Model PSdsq2J 
showed no significant change in disk structure due to the 
inclusion of a polar jet. 

In the Vt/,/2 models accretion rates were consistently 
> 0.30 M Q s" 1 and L v > 100 - 150 foe s" 1 during the 
spiral phase, with peak magnitudes occuring in the pres- 
ence of global modes. Both shellular models (full and re- 
duced Eth) had a higher average accretion rate than the 
cylindrical model, which averaged and peaked at M ~ 0.1 
and 0.4 Mq s , respectively. These models all produced ac- 
cretion rates and neutrino luminosities high enough to fuel 
LGRBs. 

In cases with rapid accretion, variations of up to half 
an order of magnitude in M and L„ occured as rapidly as 
the snapshot interval of 0.02 s. Such fluctuations have been 
observed in the spectra of LGRBs, although several factors 
in the system may combine to produce these, such as jet 
structure, environment, etc. Here, the bursts of accretion 
are due in large part to the presence of dense clumps which 
have been created by the spirals and fragmentation. 

We also note that while M and L„ show strong (pos- 
itive) correlation, no direct scaling appeared between the 
quantities. In general, the rates of the former appear to make 
the B-Z mechanism a significantly more likely scenario than 
^-annihilation for producing LGRBs, keeping in mind the 
simplifying assumptions made in calculating both. The neu- 
trino rate was probably too low much of the time, as L v 
only reached above 100 foe s _1 at brief intervals for accre- 
tion rates typically of > 0.5 Mq s _1 . 

The effects of general relativity in the system were non- 
negligible. During the evolution of the Newtonian models, 
typical (though not peak) accretion rates and neutrino fluxes 
were roughly 15 — 20% higher than in the pseudo-potential 
models. At the intermediate central object spin rates present 
in these collapsar simulations (a w 0.5), the overall differ- 
ence in potential wells was quite small (w 1%, see caption 
in Table |5). The deviations in disk temperature and struc- 
ture appear to be the determining factor in the energetics. 
This stresses the importance of including GR in as realistic 
a manner as possible in simulations, and further work on 
this task is being done. 

11.4 Nucleosynthesis 

In studying NSE-produced materials in the collapsar, it was 
found that light elements dominated the disk, and that, dur- 
ing the early evolution, the masses of 56 Ni produced in the 
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corona and beyond were too small to power a HN. Typi- 
cally the amount of 56 Ni present in the system remained 
constant after the early phases of the disk (t > 0.30 s), 
when the corona had reached an approximate state of equi- 
librium. More 56 Ni was produced in the lower v$ models, 
with the greatest amount in the low-T model, 0.10 Mq. In 
investigating areas for further production of radioactive Ni, 
the v<f, /2 models produced the strongest equatorial outflows, 
even though those disks were the smallest; PSd2Kd5 pro- 
duced the highest outward velocity. 

It was also found that, by the inclusion of a simple col- 
limated jet, rapid polar outflows were produced (here, with 
v ~ 0.1c, but this may be strongly model-dependent), which 
carried high entropy material with moderate and low elec- 
tron fraction outward. This material from above the disk and 
the coronal edges which was swept outward may be suitable 
for forming 56 Ni, and in some cases, for the r-process to 
manufacture other heavy (neutron-rich) metals. Due to the 
necessarily high velocity of any outflow, one might expect 
that, once begun, a relatively large amount of 56 Ni and Fe- 
group metals could be deposited outside the disk in a short 
period of time. 

Studying the relative locations of Fe and O did not pro- 
vide overwhelming evidence for or against the presence of 
polarised asymmetry, as has been observed and predicted 
in HNe. In general, O is prevalent in the outer corona, and 
Fe, outside that region. However, in the presence of the jet, 
material from the corona would be carried up along the po- 
lar axis, with a strong likelihood of heavy elements being 
manufactured. 

11.5 Collapsar implications 

While not spanning the lifetime of a typical LGRB, all 
models in this study ran for a large number of dynamical 
timescales, and they provide a basis for postulating con- 
tinued behaviour for long-duration phenomena. It is worth 
considering the likely behaviour of our models continuing 
after the end of our simulations, particularly in connection 
with its implications for both LGRBs and HNe. 

Regarding the former, the long-term stability and be- 
haviour of the collapsar disk must be considered. In the 
scaled models, high accretion rates occured, and the 

size of the spiral structures grew globally; in the v$ / v2 mod- 
els, the densest material (p > pdrip) tended to be accreted, 
most likely requiring reformation of the disk to permit high 
accretion rates. In all of these systems, however, the Qt val- 
ues generally remained low at late times, and a high degree 
of non-axisymmetric structure remained even in the catas- 
trophic cases, along with dense clumps of material. More- 
over, the flattening of specific angular momentum profiles 
(shown for shellular profiles in Fig. 132ft for models with de- 
creased v<f, also suggests continued instability for the disks. 
With sustained infall, the longterm disk behaviour can then 
be broadly classified into two possible scenarios: alternating 
destruction and reformation, and a continual presence with 
more 'balanced' mass loss/gain. Again, these may be evi- 
denced in and also explain some of the variety of observed 
LGRB lightcurves. While short term fluctuations would be 
due to clumps and local density differences in the flow, longer 
term ones may arise from the oscillating behaviour of the 
disk as a whole. 



The details of this behavior may strongly affect HNe as 
well, in considering the effects of time-integration on heavy 
element production. The yields from the models in this study 
must be considered as minima, as the inclusion of non-NSE 
nucleosynthetic processes will certainly lead to increased 
amounts of both 56 Ni, which powers any HN lightcurve, 
and other Fe-group products, which characterize the non- 
spherical geometry of the event. In the presence of near-jet 
outflows, even NSE yields may be enhanced as heavy metals 
are produced and ejected. Disk evolution, whether proceed- 
ing by destruction and reformation or by balanced mass gain 
and loss, will influence whether metal production is quasi- 
continual or -discrete, or whether potential Fe-group mate- 
rial is accreted rapidly, delaying any possible HN event. 

Within the scope of these trends, there is also, impor- 
tantly, a wide range of properties obtained from this related 
sample of LGRB candidates. For example, very different 
structures of the M and L v curves were obtained; equa- 
torial outflows of various velocity ranges were produced; the 
amount of heavy metals (particularly 56 Ni and Fe), synthe- 
sised in the disk and corona varied by orders of magnitude, 
even for models of similar total angular momentum. Such 
diversity exists within the observed LGRB/HN population. 
For example, a system with high accretion rates, low 56 Ni 
content in the corona and small outflows in the equatorial 
plane might produce a bright LGRB with a weak SN ex- 
plosion (due to low radioactive Ni content). Conversely, a 
system with low accretion rates and a large mass of 56 Ni 
may result in a strong HN with a weak LGRB or with a 
lower energy XRF, or possibly only a bright, polarised SN. 
We conclude that dynamically unstable collapsar disks pro- 
vide a unified and simple explanation for the large variety of 
these events. 

Among the family of potential collapsar progenitors in- 
vestigated here, the principal discriminating factor appeared 
to be the overall amount of rotation initially present. Sim- 
ulation results led to the classification of LGRB candidates 
by a simple yj = jo/iisco parametrisation with a window 
of values for systems that evolved hydrodynamically unsta- 
ble disks able to fuel a central engine. Importantly, this yj 
prescription is quite general. Given the diversity of observed 
LGRBs, the heterogeneity of the GRB-SN/HN connection 
and the probable relations with phenomena such as XRFs, 
it is likely that a variety of progenitor scenarios will require 
examination, including both binary and single-star systems, 
to fully understand LGRBs. Further considerations may re- 
fine progenitor characterisation, but we expect the trends of 
this study to be broadly applicable across a wide range of 
possible collapsars. 

Additional work should be done to further these re- 
sults, in particular by including non-NSE nucleosynthetic 
processes and a more realistic equation of state as well as 
a physically-motivated jet. These physical inputs are neces- 
sary for understanding not only LGRBs themselves but also 
their particular connection with HNe. Finally, while known 
to be essential for driving burst jets, magnetic fields were 
not directly included in this study, and their full role during 
collapse and disk evolution must be investigated further. 
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APPENDIX A: DISSIPATION OF MRI, 
B-FIELD CONDITIONS 

In c onsideration of the on set of MRI in an accretion 
disk. lBalbus fc Hawlevl (|l99ll ) estimated conditions for non- 
ideal MHD effects to damp the growth of magnetic modes 
from dimensional analysis. Specifically, the local dynamical 
timescale was compared to the diffusion timescale of modes 
due to finite resistivity and thermal conductivity. This re- 
sulted in the following criterion for damping to be negligible 
(for a z-component): 



< 1, 



(Al) 



where fix is the orbital frequency assuming Keplerian ro- 
tation; iiAz, the poloidal Alfven velocity; and x, a general 
diffusion coefficient. They stated that the minimum B-field 
strength for damping to be ignored, and therefore for the 
MRI to grow nonlinearly, was extremely small for most ac- 
cretion disks; as a result, a magnetic field, internal or ex- 
ternal, of any size would lead to the development of global 
modes in the disk on dynamical timescales. We investigate 
the same criteria in the realm of the collapsar disk structure, 
which has already been noted to be of an extreme nature. 

In considering finit e electrical res istivity, plasma theory 
leads to an expression (|Spitzerlll965P ) for the coefficient in 
Eq. CSI]): 

^ = 3.8xl0^1n(A) (A2) 



A 



7e T 3 / 2 

3 (fc B T) 3/2 _ 4.31 x 10 20 T 3/2 



2 Ze 3 {T, 



\l/2 



where Z is ionic charge; 7e ~ 1, the ratio of conductivity 
to that of an ideal Lorentz gas; and ln(A), the Coulomb 
logarithm. Substituting this coefficient into Eq. (|A1|) for a 
Z — 1 (fully ionised) gas, the general condition on the local 
poloidal B-field magnitude for dissipation due to electrical 
resistivity not to damp out the MRI is: 



B 



2.83 x 10 27 ln(A) p ( M 



(i?T) 3 / 2 



1/2 



M^ 1 



(A3) 



where M is the mass of the central BH. 

The analogous diffusion coefficient representing the 
therm al conductivity of the plasma is l|Balbus fc Hawlevl 



Xtc 



1.84 x 10" 
MA) 



T-,5/2 

-T ' cm s 



and the corresponding condition on the B-field is 



2 57.8 T 5/2 / M 
z ' tc >?> ln(A) i? 3 / 2 \M e 



1/2 



(A4) 



(A5) 



As an example, we check the requirements on the mag- 
netic field in the near-inner regions of the collapsar disk, in- 
serting reasonable values of physical quantities (R — 10 7 cm, 
p = 10 11 g cm" 3 , Y c = 0.4, T = 10 10 K, M/M = 3). The 
conditions with electrical resistivity and thermal conductiv- 
ity require that B z > 3.3 x 10 7 G and B z > 2.2 x 10 7 G, 
respectively, in order that dissipation is unable to damp 
out MRI growth. Therefore, the minimal necessary magnetic 
field strength for nonlinear growth (B ~ 10 9 G), while not 
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unphysically high, is certainly not negligible, as it is gener- 
ally for lower density disks. 



APPENDIX B: RELATIVISTIC 
PSEUDO-POTENTIALS 

The generic task of a pseudo-Newtonian potential] is to 
capture important relativistic aspects without the compu- 
tational demands of full GR. The most well-known exam- 
ple of such an approximation for a rotation-dominated flow 
is the Paczyhski-Wiita pot ential for a non-rotating BH, 
$Pw(iJ) = -M(R- 2r g ) _1 (|Paczvriski fc Wiitalll980h . and 
there are also similar potentials for motion in the Kerr met- 
ric. The accuracy of such potentials can be tested by the 
comparison of output to exact GR solutions, such as for the 
location of the innermost stable circular orbit (Risco); the 
Keplerian orbital velocity (Qk) and the epicyclic frequency 
(k), which are no longer equal, as they are in the Newtonian 
regime; and the structural and temporal properties of evolv- 
ing accretion disks, such as specific energy and dissipation. 

We note that this involves identifying the Newtonian R 
in all of these terms with the Boyer-Lindquist (B-L) coordi- 
nate, Rb-l- To quantify the resulting deviation involved in 
this, we compare the proper circumference at Rb-l ~ r- m 
to the corresponding 'Newtonian' value, 27rrt n , in the equa- 
torial plane. Examples of (a, Cb-l/Cn) are (0,1), (0.5,1.01) 
and (0.9,1.13). The ratio quickly and monotonically ap- 
proaches 1 with increasing radius, e.g. for a = 0.9, it deviates 
by only 1% at R = 3n n - 

For the cases considered here (corotating orbits and 
not extremely fast BH spin, 0.4 ^ a ^ 0.9), most pseudo- 
potentials are essentially Newtonian for radii ^ 20 — 25 r g ~ 
3 — 4 x 10 6 cm, while also reproducing the correct loca- 
tions of iiisco, which is used as the BH accretion boundary 
here. However, no single potential is able to reproduce all of 
the above-noted properties simultaneously to errors within 
even 20% in the relativistic regimes, and therefore one must 
choose to focus on certain properties over others depending 
on the system. 

We have considered four different pseudo-potentials 
for corotating motion around BHs, which we im- 
plement as acceleration terms. The asEP acceleration 
which is implemented i n th is study was derived by 
iMukhopadhvav fc Misral (|2003T ), and is given in Eq. (0 
above. As shown here, this form reproduced most exact GR 
properties within tolerable limits. 

From lArtemova. Biornsson fc Novikovl |l99(f) , the ac- 
celerations Fb and i*6 (their Eqs. (13)-(18)) are practically 
identical, and so we implemented only the former: 



Orbital Frequency (£2 K ), a B h = 

0.50 



ri 



= 1 + 

P = 



ri 



where n denotes the location of the event horizon (with ra- 
dius and spin in units of M). Figs. 4-6 of lArtemova et ail 



5 The same notations as in i|4.2l are used here. 




Figure Bl. A comparison of orbital velocities for different 
pseudo-Newtonian accelerations, including the exact GR value, 
for a=0.5. (R is measured in units of M.) 
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Figure B2. A comparison of elliptical frequencies for different 
pseudo-Newtonian accelerations, including the exact GR value, 
for a=0.5. Note: aArts and aMKep produce nearly the same fre- 
quencies for this BH spin. 



i|l996h show the use of this term in reasonably approximat- 
ing an optically thick a-disk with a low accretion rate by 
comparing optical depth (oc surface density) and tempera- 
ture to a fully relativistic case for different spins. However, a 
comparison of angular and particularly epicyclic frequencies 
in Fig. IB II and Fig. IB2I respectively, shows great disparities 
with the exact GR results for a = 0.5, and these increase 
with a. 

Similar results are pr oduced with the acceleration given 
by iMukhopadhvavl (|2002r ) , derived from exact GR solutions 
for specific energy and angular momentum. His Eq. 12, 



aMKop(-R) 



M(R 2 - 2oVR+ a 2 ) 2 
R 3 (R 3 / 2 -2^R+a 2 ) 



(B2) 



reproduces GR energy dissipation rates calculated very well 
(his Fig. 1), but does not reproduce well the rotational kine- 
matics, as seen again in Figs. IBHB2I 

In order to obtain more ac curate epicyclic frequencies, 
IMukhopadhvav fc Misral |2005t ) integrated a function ap- 
proximating the GR value of k, yielding their 'Logarithmi- 
cally Modified Potential' (LMP), which has a radial accel- 
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Figure B3. A comparison of pseudo-Newtonian accelerations, for 
a=0.5. In these inner regions, a^rts and ajyiKcp have noticeably 
larger accelerations than the other potentials, of a factor 2-3 at 
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eration: 



aLMp(-R) 



Mr m ( 9 



R 3 
M_ 



20 



■ 111 



R 



(3R- 



)2/9 



+ 



(B3) 



show that the kinematic properties of 
are fairly accu rate; however, Fig. 5 of 
iMukhopadhvav fc Misral (|2003h shows that the specific en- 
ergy of these orbits is much greater (> 30%) than the GR 
values. In contrast, the chosen potential in this study, SEP 
from the same authors, varies only slightly from the GR spe- 
cific energy (< 5%), and also from the angular frequency as 
shown in Fig. IB II above. It deviates the most from the GR 
epicyclic frequencies, up to 65% as shown in Fig. IB2I 

In the case of the collapsar system, we are studying 
global hydrodynamic stability which depends on local val- 
ues of k, as reflected by the form of Qt- The structural 
energetics are important globally, as stability is determined 
also by cooling rates, determined by temperature and den- 
sity. Furthermore, LGRBs powered by neutrino-annihilation 
require large amounts of neutrinos produced in the inner- 
most regions of the disk, where density and temperatures 
are highest. 

The relative values of the radial accelerations them- 
selves vary quite a lot, by more than a factor of 2 at the 
ftsco, as seen in Fig. IB3I The inner boundary condition 
provided at this point by a,n(R) is also key to the accretion 
rate which drives the LGRB. The GR value is calculated in 
the equatorial plane for a particle at rest in the Kerr metric, 
and is well approximated by the LMP and SEP accelera- 
tions for a — 0.5. However, it is interesting that for a = 0.9, 
aGR,(Risco) has increased much more rapidly than these, 
and is best approximated by the Art5 formulation to within 
15%, being now a factor of two greater than the LMP/SEP 
values. 

However, even at this high value of a, the SEP pro- 
duces a very accurate value of Qk throughout most of the 
inner relativistic region, unlike Art5 and MKep. The SEP 
(Art5/MKep) potential gives a peak value for k roughly a 




p log(g cm 3 ) 

Figure CI. 56 Ni abundances in NSE, as percentages of relative 
mass fraction in the relevant (p, T) phase space, with successive 
shades representing order of magnitude intervals in log(%). 



factor of 2 (3) greater than the GR one, and even the LMP 
value is ~ 50% above the GR value. 

However, over this whole range of BH spins, the SEP 
provides accurate values of specific energy (and of Qk) - Even 
though we are explicitly focused on the hydrodynamic sta- 
bility of the LGRB disks and the value of the parameter 
Qt, we have chosen a potential whose greatest weakness is 
in reproducing the kgr- However, we note that, through the 
early evolution of the disk, the minimum values of Qt, de- 
noting those regions approaching hydrodynamic instability, 
are at radii of order 5 — 10 x 10 6 cm (~ 20M), where the 
Newtonian approximation is quite good. Also, we reiterate 
the importance of the inner, relativistic regions in repre- 
senting accurate structure and energetics of the disk. We 
therefore use this SEP approximation in this study, whilst 
acknowledging its limitations. 



APPENDIX C: NSE ABUNDANCES 

All NSE calculations of A and element and isotope 
abundances used in this study have been made by F. 
Timmes using a 3000+ isotope reaction network. The 
code used is related to the publicly available NSE code 
at 



http://cococubed.asu.edu/code_pages/nse.shtml, but in- 



cluding a larger network of isotopes and more detailed 
physics, such as Coulomb corrections. Figs. ICHCllI show 
the mass fraction abundances of the NSE material analysed 
in this work, in the relevant (p, T) phase space. 

This paper has been typeset from a T^X/ F/IFX file prepared 
by the author. 
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Figure C6. Ca abundances (and see Fig, IClt . 
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Figure C8. O abundances (and see Fig. ICl|l . 
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Figure C7. Si abundances (and see Fig. ICl| l. 
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Figure C9. 4 He abundances (and see Fig. IClj . 
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Figure CIO. Proton abundances (and see Fig, ICll l. 




Figure Cll. Neutron abundances (and see Fig. ICll l. 



